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SUMMARY 


Accounting for the statistical geometric and material variability of structures in analysis has been 
a topic of considerable research for the last 30 years. The determination of quantifiable measures of 
statistical probability of a desired response variable, such as natural frequency, maximum displacement, 
or stress, to replace experience-based “safety factors” has been a primary goal of these studies. There 
are, however, several problems associated with their satisfactory application to realistic structures, such 
as bladed disks in turbomachinery. These include the accurate definition of the input random variables 
(rv’s), the large size of the finite element models frequently used to simulate these structures, which 
makes even a single deterministic analysis expensive, and accurate generation of the cumulative 
distribution function (CDF) necessary to obtain the probability of the desired response variables. 

The research presented here applies a methodology called probabilistic dynamic synthesis (PDS) 
to solve these problems. The PDS method uses dynamic characteristics of substructures measured from 
modal test as the input rv’s, rather than “primitive” rv’s such as material or geometric uncertainties. 
These dynamic characteristics, which are the free-free eigenvalues, eigenvectors, and residual flexibility 
(RF), are readily measured and for many substructures, a reasonable sample set of these measurements 
can be obtained. The statistics for these rv’s accurately account for the entire random character of the 
substructure. Using the RF method of component mode synthesis, these dynamic characteristics are used 
to generate reduced-size sample models of the substructures, which are then coupled to form system 
models. These sample models are used to obtain the CDF of the response variable by either applying 
Monte Carlo simulation or by generating datapoints for use in the response surface reliability method, 
which can perform the probabilistic analysis with an order of magnitude less computational effort. Both 
free- and forced-response analyses have been performed, and the results indicate that, while there is 
considerable room for improvement, the method produces usable and more representative solutions for 
the design of realistic structures with a substantial savings in computer time. 
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TECHNICAL MEMORANDUM 


DEVELOPMENT OF A PROBABILISTIC DYNAMIC SYNTHESIS METHOD 
FOR THE ANALYSIS OF NONDETERMINISTIC STRUCTURES 


1. INTRODUCTION 


1.1 Motivation 

The primary assumption for the basic analysis of structures has always been knowledge of the 
material properties and of the geometric layout. Once an initial analysis is performed, uncertainties in 
the above inputs, as well as uncertainties in the analysis method used, computer accuracy, and many 
other factors are addressed by the use of somewhat arbitrary experience-based “safety” factors. Some 
attempt has been made over the last several decades to account for the variation in material properties by 
using “3— a” curves from material handbooks, but in general, the statistical variability of structures has 
not been quantitatively accounted for in structural analysis. 

This problem has been well recognized and therefore a significant amount of research has been 
performed over the last 30 years to develop methods of dealing with the statistical variability of struc- 
tures. The goal of these methods is that given statistically characterized input quantities, a desired output 
response quantity can also be characterized statistically. Examples of the varied input quantities, which 
can be defined as random variables (rv’s), are Young’s modulus, density, thickness, length, and width. 
Output rv’s could be, for example, stress, load, or natural frequency. There are various levels of statisti- 
cal characterization for these rv’s; the most complete is the cumulative distribution function (CDF), 
which defines the probability of occurrence of a specific value of an rv over it's entire range of 
possibility. 

For very small problems of one or two input rv’s, an exact solution for the response rv can be 
obtained. The first step is to obtain the joint probability distribution function, which is the multidimen- 
sional integral defined by the probability distribution function (PDF) of the input rv’s. This surface (or 
curve for one rv) must then be integrated from negative infinity to the domain defined by the function 
relating the input to response variables. For multiple degree-of-freedom (dof) structures, however, this 
integration often becomes intractable. For these multi-dof structures, the most straightforward method of 
obtaining these results is to use the Monte Carlo (MC) method. In this method, a complete sample set of 
each input rv fitting a desired distribution is created. A variety of methods are available for generating 
these sample sets. 1 Once a set for each input rv is created, an analogous solution sample set is calculated 
using a numerical method capable of analyzing large structures, such as the finite element method 
(FEM) or boundary element method. The statistical parameters of the desired response value is obtained 
from this solution set and a distribution type can be fit to it using commercially available software. 2 The 
drawback to this method is the large number of runs required to obtain a response CDF and the even 



larger number required to obtain an accurate PDF. Although the capability for performing these calcula- 
tions for realistic multi-dof structures is improving dramatically and will be discussed in this disserta- 
tion, the process is still extremely computer-intensive and expensive, especially if multiple runs are 
required for parametric design studies. 

The thrust, therefore, of the published research has been to develop approximate methods, 
frequently called reliability methods, to obtain response statistics of complex structural components 
using an order of magnitude less calculations than MC techniques. Initial methods assumed Gaussian 
input rv’s and obtained only the mean and standard deviation of the response rv. Further development 
has resulted in advanced, commercially available computer codes, such as Southwest Research 
Institute’s NESSUS code , 3 that allow input rv’s of any distribution type and produce entire distributions 
of response variables along with accurate statistical parameters from finite element structural models. 
These codes incorporate a number of reliability methods that reduce the number of runs to as few as the 
number of input rv’s plus one, down from the several thousand runs required for MC. 

There are still many limitations to the direct application of the above codes to complex structural 
components. One main problem is determination of the statistical characteristics. of the input rv’s. A 
great deal of data are available characterizing material property variability, but accurately defining 
geometric variation in a structure is extremely difficult. In addition, because the structure is typically 
represented by spatial coordinates of finite element grids, using the location of each of the grids as an rv 
would make the probabilistic problem intractable. The method of random fields can be used to describe 
the spatial variation using correlation techniques, but obtaining accurate correlation for a manufactured 
part is frequently not possible. 

Another problem is that many of these structural components require finite element models of so 
many dof ’s that substructuring reduction techniques are required to decrease the size of the model just to 
perform a single, deterministic analysis. These techniques separate the model into groups of elements 
called substructures; frequently these substructure models are composed of similar groups, such as 
turbine blades on an engine turbine bladed-disk, or are required because different companies build the 
substructures, such as the booster, external tank, and orbiter of the space shuttle system. The mass and 
stiffness matrices of these substructures are then mathematically reduced to a smaller set; one group of 
reduction techniques called component mode synthesis (CMS) methods represent the substructures by a 
set of interface dof ’s and generalized modal coordinates. The combination of CMS methods and reliabil- 
ity probabilistic techniques is one of the topics explored in this thesis. 

One possible solution to the above problems is the direct use of statistics generated by modal 
testing of a population of the substructures. One example of such a population is turbine blades for a 
specific stage in a turbopump. For many engine development programs, the experimental determination 

a hundred blades each, the population can therefore number in the thousands. The measured values from 
this testing, which include the structural eigenvalues, eigenvectors, and residual flexibility of the bound- 
aries, would be established as the “dynamic” rv’s of the system, rather than “primitive” rv’s like variable 
thickness or length. These dynamic rv’s can be used directly in the residual flexibilty (RF) method of 
CMS to create the substructure mass and stiffness matrices. 
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The objective of the research presented here is to apply the RF CMS method to the probabilistic 
analysis of complex structural systems using this set of statistically accurate “dynamic” rv’s as input. 

The method developed is called probabilistic dynamic synthesis (PDS). Application of this new method 
would enable accurate probabilistic analyses to be performed on a variety of realistic, complex structural 
components that presently are not performed due to enormous computer requirements and limited result 
accuracy. One type of structure particularly applicable to the PDS method is the mistuned bladed disk; 
the size of blade finite element models is so large that realistically modeling even a single sample of 
such a disk, with every blade slightly different, has not been attempted even with CMS techniques. Only 
recently have studies using reduced order methods on bladed disks been initiated. 4 The use of PDS 
allows a complete probabilistic analysis to be performed of the mistuned bladed disk to generate both 
free- and forced-response solutions. 


1.2 Scope of Dissertation 

Because one of the major objectives of this research is to combine several different disciplines 
together to form a new analytical tool, extensive theoretical and literature background will be presented 
on each of these aspects. The first topic to be examined is the CMS procedure, which is presented in 
chapter 2. The emphasis in this chapter is the comparison between the CB and residual flexibility meth- 
ods of CMS and why residual flexibility is more applicable to this research. 

Chapter 3 is an extensive study of probabilistic analysis as applied to structures. This study is 
divided into six parts: (1) theoretical basis and exact methods, (2) MC approximate method, (3) reliabil- 
ity approximate methods, (4) comparison of reliability methods, (5) perturbation methods, and 
(6) comparison of methods and their areas of applicability. 

Chapter 4 presents the application of PDS to a simple test case consisting of two substructures 
represented by masses and springs. Some of the remaining aspects of the theoretical background of PDS 
will be completed as part of the examination, including examining the effects of correlation of the input 
rv’s. 


Chapter 5 examines the application of the PDS method to a realistic three substructure system 
modeled using the finite element method. Both free- and forced-response solutions are presented, along 
with examinations of many aspects of the problem that are only encountered because of the realistic 
aspect of the models. In addition, an MC analysis of the system is described; this type of analysis on a 
structure with this level of complexity has not previously been reported in the literature. 

Chapter 6 applies the PDS method to realistically modeled tuned and mistuned bladed disks. A 
literature survey of bladed disk research will be presented, and the results from chapter 3 will be used to 
show the advantages and disadvantages of using PDS for this application. Improvements in the method- 
ology resulting from the work presented in chapters 4 and 5 are applied, and both free- and forced- 
response solutions are produced. 

Chapter 7 discusses the applicability of the method to general problems in structural analysis. A 
detailed chart is presented that lists the methods used in this field of research and summarizes the advan- 
tages and disadvantages of each. Some specific directions for future research are examined, and final 
conclusions of the research are presented. 
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2. THEORETICAL BACKGROUND AND LITERATURE SURVEY— COMPONENT 

MODE SYNTHESIS 


CMS is a very well developed topic in structural dynamics. An extensive literature survey will 
not be presented here; rather, the applicability of the various methods to PDS will be discussed and the 
theoretical background behind those procedures will be presented. CMS is extensively used in the 
automotive and aerospace fields, and was initially discussed by Hurty in 1965. 5 In 1968, Craig and 
Bampton developed the concept into a readily usable methodology and it became widely used by indus- 
try. 6 The CB method uses fixed-interface modes and “constraint” modes to represent each substructure. 
An alternative method, initially investigated by MacNeil, uses free-interface modes and residual flexibil- 
ity of the interface to represent the substructure instead. 7 In 1984, Martinez altered the method into a 
representation similar to the CB method for ready use by industry. 8 

The CB method was initially evaluated for its applicability to the PDS method. In this method, a 
structure is divided into substructures a,b,...,n, which are denoted with superscripts sa through sn. The 
physical displacement vectors of each substructure, which have either a subscript i denoting internal dof 
or a subscript b denoting boundary dof, can be written as 



Similarly, the substructure mass and stiffness matrices can be written in partitioned form as 



M ib ' 

sa 

'Kji 

K ib ' 

sa 

"Mjj 

M ib ' 

sn 

r Kii 

X 

_M bi 

M bb_ 

9 

1 

2"! 

K b b. 

9 • • ■ 9 

M b i 

M bb. 

9 

i 

cr 

K b b_ 


The CB transformation between the physical coordinates and a reduced vector of generalized and physi- 
cal coordinates is 


x i 

x b 



= [cb] sa 



The partitioned transformation matrix [cb] sa is defined as 


[cb] sa 


O n ¥ 
0 I 


(3) 


(4) 


4 



where the “constrained transformation modes” matrix is defined as 


m = -[K ii ] sa ' 1 [K ib ] sa . (5) 

and O n is the truncated matrix of N cantilevered modes of the substructure (the CB method uses 
clamped boundary conditions, as mentioned previously). The vector of generalized coordinates in 
equation (3) is composed of q, a truncated set of generalized coordinates that describe, approximately, 
the internal displacement field, and x b , the physical coordinates of the substructure boundary that are 
retained throughout the analysis. The size of q is equal to the number of modes retained for the analysis. 
The number of retained modes depends on a variety of factors including the desired accuracy, the fre- 
quency range over which the model is intended for use, the frequency range of the input excitation, etc. 
Generally the number of retained modes is much smaller than the number of physical dof ’s of the sys- 
tem and is the source of the reduction in the size of the problem. 

The transformation between each substructure’s reduced vector and the total structure’s reduced 
vector is then calculated. For substructure a, this transformation will be 


where 


q 


sa 





[T] 


sa 


... 0 

... 0 


0 ••• 
I ••• 


o' 

0 • 


( 6 ) 


(7) 


Using equation (3), the transformation between the system-reduced vector and the substructure physical 
vector is therefore 
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where 


X; 


i sa 


l x bj 


= [«] 


isa 


,sa 


sn 


v sa 

x b 


I Y Sn I 

l x b J 


( 8 ) 


[a] sa = [cb] sa [T] 


i sa 


( 9 ) 


A conservation of energy approach is used to synthesize the system mass and stiffness matrices 
for the system undamped equations of motion. A complete derivation is given in the textbook by Craig. 9 
Essentially, the kinetic and potential energy of substructures are defined as 


[K.E.] sa = 



[P.E.] sa 



( 10 ) 


The energy terms for each substructure are then added together to form the total kinetic and potential 
energies of the system, which can be expressed in terms of the system coordinates using equation (8): 
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'f 

T 

q sa 


q“ 

T 


[K.E.] sys = -j* 

q sn 

x S b a 

> [a] saT [M] sa [a] sa < 

q sn 

x S b a 


q sn 

^.sa 

x b 

> [a] snT [M] sn [a] sn < 

q“, 

x sa 

x b 






Y Sn 

l x b J 


Y SI1 

l x b J 


f 

• Qfi 

T 


q 


q sa 

/. sn 


• sn 

q 

• sa 

> [M] sys < 

1 

* sa 

x b 


Xb 

• sn 


• sn 

l x b J 


l x b J 



q sa 

T 

q sa 


q sa 

T 

q sa 

[P.E.] sys = 

q Sn 

x S b a 

• [a] saT [K] sa [a] sa < 

q sn 

sa 

x b 

> + ...+A< 

q sn 

v sa 

x b 

► [a] snT [K] sn [«] sn < 

q sn 

x S b a> 


x b\ 


sn 

l x b J 


Y sn 

l x b J 


x s b n . 



T 

r ^ 

e a 

q a 


q sa 

„sn 


^sn 

q 

sa 

> [K] sys < 

q 

sa 

x b 


x b 

„sn 


sn 

l x b J 


L x b J 


7 



where 


[M] sys = [«f aT [M] sa [a] sa + • • • + [a] snT [M] sn [a] sn 
[K] sys = [«] saT [K] s >] sa + • • • + [a] s " T [K] sn [a] sn 


These matrices can now be used to form the system equation of motion, 


[M] 




q sa 

~sn 

> [K] sys^ 

> 

^sn 

q 


q 

.V 


. x b. 


( 12 ) 

(13) 


(14) 


from which the system eigenvalues and eigenvectors can be obtained. 

The CB procedure is a straightforward and time-tested method of substructure synthesis. How- 
ever, upon closer examination, it was determined that the CB method relies upon data that are not easily 
obtained through experimental tests. This can be seen if one applies the CB transformation as defined in 
equation (4) to a single substructure: 10 


M sa l cb = 

cb sa ] n 

[ [M sa l[cb sa l = 

I 0 N (M ii T / + M ib ) 

J 

* 

JL J 

!_sym T' T (M ii T' + M ib ) + M ib T' + M bb _ 


(15) 



cb. 




K 


bb 


0 

K bi Ku‘K bi 


(16) 


The A n partition is the diagonal of N retained substructure eigenvalues (squares of the natural frequen- 
cies) and the other variables and subscripts refer to equations (2), (4), and (5). It is clear that significant 
portions of the matrices are derived directly from the original mass and stiffness matrices and cannot be 
obtained from test (unlike A N and <I> N ). Recall that a major goal of this work is to incorporate random- 
ness into the model using information that can be measured solely from test. This information is not 
fully accounted for in the above transformed mass and stiffness matrices. However, since the method is 
fairly simple to use, it does suggest that it would be the CMS method of choice for a purely analytical 
PDS methodology, or possibly a hybrid that would use some test information and other analytical infor- 
mation. This might be useful for parametric mistuning studies, for instance. 

In order to circumvent the limitations of the CB method, the RF method was investigated. One of 
the fundamental concepts used to reduce the dof’s in CMS is the truncation of substructure modes at 
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some level, thus assuming that the higher modes will not have a major effect on the combined system 
modes. The RF method partially accounts for the higher modes by determining their flexibility. A side 
benefit is that all the elements of the system stiffness matrix can be obtained from test and that the mass 
matrix can be closely approximated by a unity matrix in the nonboundary partition (equal in size to the 
number of kept modes N). Since all the substructure information can be obtained from test, probabilistic 
data can be completely incorporated into the system matrices to obtain the system modes. 

The derivation of the free response of substructures using the RF method is presented below. The 
presentation, which is similar to that of Admire and Tinker, 11 begins with the statics equation for the 
structure: 

[K]{x} = {F} (17) 

where [K] is the stiffness matrix and { x } and {F} are the nodal displacement and force vectors, respec- 
tively. Solving for { x } , the flexibility formulation is obtained: 

{ x } = [G]{ F } where [G] = [K] _1 • (18) 

At this point, a modal analysis is performed and the lower N modes are retained for future use. The 
vector {x} can be partially represented by a transformation to the generalized coordinates {q N } using 
these N retained modes of the substructure 

i x partial}=[ <1>N ]{q N } • (19) 

Substitution of equation (19) into equation (17) and premultiplipation by [<f> N ] T yields 

® N ] T [K N ][<I> N ]|q N ) = [4. N ] T {F) . (20) 

If the retained modes O n are mass-normalized, 12 this is equivalent to 

[A N ]{q N } = [o N ] T {F| . (21) 

Now, using equations (19) and (21), the following result is obtained for {x} partia j: 

! Xpmtial 1 = [<t> N ][ a N ] [<1> N ] T ( F I = [a retained ] i F ( . (22) 

The total residual flexibility [G] can therefore be split into two parts, the first resulting from the retained 
modes (shown in equation (22) and the second from the deleted modes, which is defined as the residual 
flexibility [G residual ]; 
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{ X } — ([G re tained] [^residual 


(23) 


and [G residual ] can be solved for as 


[Gresidual] = [°] - [® N ][a N ] [® N f • 


(24) 


Therefore, the physical displacement set can be completely represented by two terms: 


x = 


® N ]{q N } + [G re sid„ai]{F} . 


(25) 


As with the CB method, the substructure is now divided into internal dof ’s and boundary dof ’s: 



(26) 


For a free-response analysis of a structure, the only forces that must be accounted for are the internal 
forces that act at the interfaces of the substructures. Therefore, the equation of motion for the free- 
response analysis of the substructure is 


Mu M ib 
M bb 

Since the applied forces at the internal dof ’s of each substructure are zero, only the boundary columns of 
the residual flexibility matrix are needed, which is beneficial for coupling with the other substructures. 
Equation (25) therefore reduces to 


• + 


l x b. 


'Ky K ib 
_ K bi K bb 



(27) 


x i 

x b 


<t>b 


^res— ib 
^res— bb 



Using the lower partition to solve for { F b } yields 

{FJ = G;‘ _ bb x b - G;‘_ bb O b N q N . 
Substituting equation (29) into the top partition of (28) yields 


x i - q N + G res _i b ^G r J s _ bb x b - G r e S _ bb O b q N j 
= — G res _ ib G res _ bb <I> b jq + G res _j b G res _ bb x b . 


(28) 


(29) 


(30) 
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Substituting this result into equation (28) yields the RF transformation matrix [T RF ] 


x i 


l x b. 




N 


G res-ib G res-bb®b 

0 


G res— ib G res-bb 


qf 1 1 _ t t rf 


x b 


7 {qil 


(31) 


Applying this transformation to the undamped equation of motion (24) to obtain the equation in the 
reduced coordinates |qj} yields 


M 


RF 


x b 





(32) 


where 



and 


K 




* N i N F - 1 N 
A +<P b Cj res _ bb <P b 


sym 




-1 

res-bb 


-1 

res-bb 


(33) 


(34) 


and [M] and [K] are the partitioned mass and stiffness matrices from equation (27). The symbol O b N 
refers to the partition of the modal matrix consisting of N retained modes on the boundary of the sub- 
structure. These substructure matrices are then coupled together to form the system mass and stiffness 
matrices in an identical manner as used in the CB method. 


To facilitate ease of implementation, the approximation shown in equation (34) is used in this 
method. The “residual mass” terms on the off-diagonals have been neglected, which makes the creation 
of the substructure mass matrix considerably simpler; this assumption is discussed in detail by Admire 
and Tinker 13 and is shown to be acceptable. In addition, the boundary partition terms on the diagonal in 
the mass matrix that are equal to zero are instead set to very small numbers in the analysis to provide 
numerical stability for eigensolution routines. 

The applicability of the RF method to probabilistic dynamic synthesis can now be seen by 
examining the substructure matrices. All of the elements of the substructure stiffness matrix can be 
obtained directly from modal tests. The free-free modes and mode shapes are generally readily available; 
if it is difficult to obtain them (for very small structures like turbine blades, for instance), recent studies 
by Han, Chen and others have identified a method for obtaining the free-free modes using constrained 
modal test data. 14,15 The residual flexibility partitions are obtained by measuring the slope of the stiff- 
ness line for drive point frequency response function measurements. The determination of this value has 
been the subject of considerable research in the last few years. 16 
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An assumption in this application of the RF CMS method is that the modal test measurements 
are accurate. There has only been limited success in synthesizing mass and stiffness matrices from test, 
and one reason may be that there are inherent errors in the measurements themselves. The fact that an 
entire sample set of these measurements is required for use in the PDS method, as will be discussed in 
chapter 4, may mitigate this error. This subject will be discussed further in areas for future research in 
chapter 7. 
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3. THEORETICAL BACKGROUND AND LITERATURE SURVEY— PROBABILISTIC 

ANALYSIS OF ENGINEERING QUANTITIES 


For most engineering problems, a desired response variable is a function of some given input 
variables. As discussed previously, if the inputs can be expressed as rv’s, then the response will possess 
a probability density function, cumulative distribution function, the statistical moments of mean and 
standard deviation, and possibly other characterizing parameters. This chapter includes a brief basic 
theoretical background of the subject including “exact” methods for obtaining the statistics of the re- 
sponse variable and discussions of various approximate methods for this determination. These approxi- 
mate methods include the MC method; reliability methods, which have grown out of civil engineering 
applications of strength versus load; and the perturbation method, which has been used extensively in 
research, particularly in the vibration analysis of bladed disks in turbomachinery. Two detailed examples 
are presented to illustrate these various techniques. Finally, a detailed comparison of the reliability and 
perturbation methods will be presented. A thorough understanding of this information is necessary to 
determine which method would be best applicable for use in the PDS technique. 

3.1 Exact Methods 

Given two continuous rv’s X and Y, the joint probability density function f XY is defined such 
that the following integral can be evaluated to determine the probability that X will take on a value 
between two limits “a” and “b” and that Y will take on a value between two limits “c” and “d:” 


db 

P(a < X < b,c < Y < d) = J|f XY (x,y)dxdy ? ( 35 ) 

c a 

where x and y are the variables of integration. The probability that both X is below the variable x and Y 
is below the variable y is defined as F XY , the joint CDF of X and Y. The integration of f XY from negative 
infinity to x and y will result in F XY as a function of x and y: 


X 


F XY (x,y) = P(X<x,Y<y)= J 

— oo 


y 


J f xY( x >y) dxd y 


(36) 


Like the functions of one variable, the joint CDF varies between 0 and 1 , where 1 is the total volume 
under the joint PDF f XY . 17 

In most engineering applications, a desired response variable is some function of the rv’s. This 
response variable can also be described by either its PDF or CDF. For a single input rv where g is a 
function H(x), the probability density function of the continuous rv G can be shown to be 18 


13 



fo(g) = f x<X> 


(37) 


dx 
dg 

For two input rv’s an exact methodology has been formulated called the auxiliary variable technique. 19 
It is beyond the scope of this thesis to describe the technique in detail, but the final expression for the 
probability density function of g requires the creation of an auxiliary dependent variable so that the 
number of dependent variables equals the number of input rv’s. This allows the following relationship to 
be expressed: 


f (c . n f xv( x »y) 

t GlG2(gl’g2j- jjj , (38) 

where IJI denotes the determinant of the Jacobian matrix, 

jgl dgi 
dx dy 
dg2 <?g 2 
dx dy 

The PDF of gj can therefore be isolated by integrating out g 2 over a region R bounded by the solution of 
g 2 in terms of x and y, which are, in turn, in terms of g,. This results in 

f Gl(gl) = J XY |j| ’ y -dg2 • (39) 

This can only be solved analytically for a few simple functions g=H(x,y) combined with simple PDF’s 
for x and y, such as a uniform distribution. Any other combination requires numerical solutions, even for 
just two rv’s. 

The CDF of the response variable can be obtained directly by integrating the volume under the 
probability density function bounded by the curvilinear plane defining a specific value of the function 
g=H(x,y): 



F g (g) = P(G < g) = |J f X Y ( x * y)dxdy . (40) 

R 

This volume is shown in figure 1 . To perform the calculation, y must be expressed in terms of x and then 
x expressed in terms of g to create the region R forming the integration bounds. Again, for any but the 
most simple functions combined with uniform distributions, the exact determination of the integral even 
for just two input rv’s is not possible; defining the domains for a nonlinear function is very difficult and 
typical bivariate distributions like the binormal distribution are not readily integrable. 
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Figure 1 . Joint probability volume of two normal rv’s. 


It is clear from the above discussion that these methods are infeasible for determining the CDF of 
an rv that is a function of many input rv’s. For applications to structural systems, in which the entire 
geometry as well as material characteristics can be random, this limitation makes the integration intrac- 
table. However, the theoretical development above has been used to develop more practical methods 
which are discussed next. 


3.2 Approximate Methods — Monte Carlo 

The MC method is an approximate numerical method for calculating the PDF and CDF of a 
function of rv’s. It is frequently used as a baseline against which other approximate methods are mea- 
sured. The method essentially takes a sample from each of the input rv’s and uses them to generate a 
solution for the response variable of interest. A large number of these solutions are produced and a PDF 
and CDF are calculated from these results. This “brute-force” method is very straightforward, but it is 
extremely computationally intensive since a large number of solutions are required for accurate results. 

The first step in the MC method is to create a set or ensemble of input variables having the same 
statistics as the input rv’s of the problem in question. A digital computer is used to generate a large 
number of uniformly distributed “pseudo-random” numbers (the vector is not completely random be- 
cause it is repeatable for the same given “seed” value). An algorithm to generate an analogous vector of 
random numbers that fall into a desired statistical distribution type is then applied. For example, the 
“Box-Muller Algorithm” is used to obtain a vector of two statistically independent normally distributed 
rv’s from a uniform set. 20 If 
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{ X }~U(0, 1 ) 


(41) 


denotes a vector X that has a uniform distribution between 0 and 1 , then 


y i = (-2 In Xj )- 5 cos (27tx i+1 ) 

(42) 

y 2 = (-2 In Xj )- 5 sin (27tx i+1 ) , 

(43) 


with 

y { and y 2 ~ N(0,1) , 

denoting that y^ and y 2 are two independent rv’s with standard normal (or Gaussian) distributions (mean 
0 and standard deviation 1). These can be converted to nonstandard normal distributions z with given 
mean |i z and standard deviation ct z using 


z = (i z + ct z y . (44) 

Once the desired distributions are created, the function of interest (e.g., t=t(zj,z 2 ) ) is simply 
calculated for each set of input rv’s z 1 and z 2 and a histogram of t is created to obtain its PDF and CDF. 
The size of the sample required depends on the CDF value of interest, which is the probability that t will 
be less than some given constant value T, denoted P(t < T). The sample size is also a function of the 
required confidence interval C, which is a measure of the percentage of the time that a sample of size N 
will yield a CDF value at P(t < T) that is accurate. One estimate of this value is 


N > 


-ln(l-C) 
P(t < T) ' 


(45) 


This necessitates a large number of samples required to obtain accurate CDF’s at small probability 
levels. 21 


3.3 Approximate Methods — Reliability 

The PDS method described in this dissertation employs the reliability method approach for 
determining the statistical structural response characteristics of interest. This methodology grew out of 
approximate techniques developed to obtain the moments of the response variable g directly from the 
moments of the input rv’s. 22 For clarity, the theory will be first presented for a response variable g that is 
a function of two input rv’s, x and y. The first two moments of each rv are defined as the means |i x and 
|i y and the standard deviations ct x and a y . These response moments can be used to determine the ap- 
proximate probability density function of g, and by themselves are useful for design and analysis. The 
first step is to expand g as a Taylor’s series about the means of the input rv’s. 
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( \ , v ^§(Mx ’ My ) , v ^§(Mx ’ My ) . N 

g(x,y)= g(M x ,My) + ^p 2L -( x “Mx) + ^^-(y^My) 


3 2 


^ S(Mx’My), v2 ^ S(Mx’My). .2 o ^8 / \/ \ 

< x ~fe> + ^2 (y-Z'y) - | - 2 ^t(’ 1 -fc)(y-M y ) 


dxdy 


( 46 ) 


Now, the mean of g can be obtained by taking the expectation of both sides 


iUg = g(^, % ) + i( ^f y) [ E (x 2 )-^] 


l g(^y)r c . ,2, „2 




[E(y 2 )-jUy 


+2 ^ E l (x " ^>( y " ^y)] ) + • • • 


(47) 


and the standard deviation of g is 


y) 


[E(xV^] + p|^j [E(y^)-^] 


+ products of cross partials + higher order terms + . . . . 


(48) 


These expansions obviously have to be truncated at some point for tractability. Wirshing performed 
detailed studies in 1975 of the errors in the mean and standard deviation resulting from using only first 
order, second order, and other selected terms from the series for selected linear and nonlinear functions 
of the rv’s. 23 As can be seen from equation (47) and (48), the error in the mean is directly related to the 
degree of nonlinearity, indicated by the second partial term, and the size of the standard deviations of the 
input rv’s, as indicated by the other term in the product. A relative measure of the size of the standard 
deviation is the coefficient of variation (cov), which is defined as the standard deviation divided by the 
mean. For some simple nonlinear functions of two rv’s, Wirshing found that the error can be up to 
20 percent using only the first order terms for high levels of the coefficients of variation such as 0.2. It is 
therefore clear that first order methods can result in significant errors for evaluating the mean and 
standard deviation in the case of nonlinear functions of rv’s. 

“Reliability” methods were developed by initially adapting the terminology discussed above to 
the analysis of structures. The goal of reliability analysis is to determine an expression for the probabil- 
ity of failure of a structure, which is one minus the probability of reliability. The exact value is obtained 
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by calculating the volume of the joint probability function of the vector of multiple input rv’s X that 
influence the structural characteristics over the failure region Q : 

P f = | £J f(X)dQ , (49) 

where P f is the probability of failure and f(X) is the joint probability distribution function for X. Due to 
the difficulty of evaluating this integral, the following approximate procedures have been developed to 
determine P f . For these methods, the function of input rv’s g is defined as the “limit state function” g(X): 

g = strength - load . (50) 

For this definition, the probability of failure is the probability that g < 0. Cornell, 24 using the reliability 
approach, developed the “first order, second moment” method by truncating the higher order terms from 
the Taylor series derived in the previous section: 


N 3 

g(X) = g(/bc) + X J~( x i -A*i) . 

i=l U 1 

This resulted in the following simple approximations for p and o of g: 

= ; 


(51) 


(52) 


r dg ^ 2 

ydx iy 


<V 


(53) 


Hasofer and Lind 25 further refined and expanded this method, now entitled the first order reli- 
ability method (FORM). They redefined g as 

g(X) = Y(X) - y , (54) 

where y is a specific value and the performance function (or response surface) Y(X) is a function of the 
rv’s. The plot of g(X)=0 will be a line on the X plane, dividing it into two parts, g < 0 (Y < y) and g > 0 
(Y > y). The probability that the function Y not exceed the value y is the probability that g < 0. For 
example, if Y(X) is a specific eigenvalue of interest, and we wish to know the probability that Y(X) is 
less than a particular value y, then that probability is equivalent to P(g < 0). There will therefore be a 
limit state function g for every value of y possible over some range of interest. If the distributions of the 
rv’s are normal and the resulting distribution of g is normal, then the following steps can be performed 
to obtain the probability that g < 0 (Y < y). First, a transformation of g to standard normal coordinates Z 
is made as follows: 
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(55) 


Z = 


g (X)- Al- 


and if P, which Cornell had termed the safety index, is defined as 


then 




_ 


P(g(X)< 0) = p(Z<-P) = O(-p) , 


(56) 

(57) 


where ® is the CDF of the standard normal distribution function found in handbooks. Since negative 
values for p are not tabulated, the relationship 

O(-P) = 1 - O(P) (58) 

is used instead to calculate this probability. The complete CDF for Y(X) is formed by finding the P’s for 
all the specific values of y in the probability range of interest. One problem with this method is that g 
can be formulated differently for some cases, thereby yielding different probabilities. To provide invari- 
ance with respect to the formulation of g, Hasover and Lind introduced an initial reduction of each of the 
primitive normal rv’s Xj into standard normal rv’s Uj using 


x i“/* Xi 

• ( 59 ) 

A i 

This transformation allows the new limit states g(U) = Y(U) - y = 0 to be plotted in standard normal 
joint probability space for every possible y. An essential point of the reliability methods becomes evident 
at this point. For a linear limit state function of normally distributed rv’s, the multidimensional bell- 
shaped volume under the joint probability density function bounded by the planar limit state function is 
exactly equivalent to a one-dimensional probability density function up to the point P (see fig. 2). The 
joint probability of g < 0 will be the volume under the joint PDF over the area where g < 0. p can be 
shown to be the shortest distance from the origin to the line g(U) = 0; the point on the line g=0 having 
the shortest distance (P) to the origin is termed the most probable point (MPP) because that point will 
have the highest probability of occurrence of any point along the line g(U)=0. The MPP is also called the 
design point X*, because it is the vector of input rv’s that will generate the probability level obtained 
from <b(-P). For a nonlinear limit state, as shown in figure 1, one source of error for this probability 
value is the degree to which the volume bounded by the curved plane is different than that bounded by 
the linear plane. Because the volume rapidly drops off away from the MPP, this error is usually small. 
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Figure 2. Reduction of two-dimensional joint PDF to one-dimension. 


The derivation above is used in the case when the origin is in the “success” region (i.e., 
g(U=0)> 0) and will yield probability levels < 0.5. If the origin is in the “failure” region (g(U=0) < 0), 
then (3 is defined to be the negative of the distance from the origin to the MPP. 26 Applying the Gaussian 
CDF function for this value will therefore yield probabilities > 0.5; if the distance 8 is defined such that 

(8 - -(3) > 0, then 


cp(-(3 ) = 0(8) > 0.5 . (60) 

An example of a closed-form limit state function is shown below to illustrate this method. Con- 
sider a cantilever beam of length L = 1 and with statistically characterized input rv’s load P and plastic 
moment strength M (see fig. 3). P is normally distributed with mean of 12 and standard deviation of 1, 
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and M is normally distributed with mean of 40 and standard deviation of 3. Failure can be defined as 
M<P*L, so a limit state function in M and P is formulated as 

g(M,P) = M-PL , (61) 

where g<0 corresponds to failure. The primitive rv’s M and P are now converted to standard normal 
variables m and p. 


(M - 40) n 10 

m = p = P-12 . 

3 


(62) 


In terms of m and p, the limit state function is 


g = (3m + 40) - (p + 1 2) 


(63) 


Setting g=0 results in an equation for m as a function of p, which can be plotted in the reliability space 
as shown in figure 3. The design point, or MPP, is located at the shortest distance from the limit state g 
to the origin. The distance (3 can then be calculated using the distance equation from the origin to a line 
and substituted into equation (58) to obtain the probability of failure for this system. 



Figure 3. Reliability space for two standard, normal rv’s. 
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In the general case, the FORM contains several sources of error for the prediction of P f , espe- 
cially if the actual response function is nonlinear. First, it assumes that the original input distributions are 
Gaussian, although this assumption is addressed. 

later in this section. Second, it approximates the slice between the failure and safe regions of the volume 
under the joint probability surface by a straight plane, as opposed to the curved surface that would arise 
from a nonlinear g function. Finally, one of the main sources of error is the inaccuracy in the calculated 
location of the MPP’s. Errors in the location of an MPP can arise due to a number of reasons, some of 
which are identified in examples described in section 3.4. 

The reliability method was expanded by Rackwitz in 1976 to multidimensional problems for 
which the limit state curve g=0 is an explicit, nonlinear function of the rv’s. 27 This, of course, makes the 
determination of p much more difficult. The method makes use of an iterative algorithm to find (3. 
Initially, an assumed value of the design point X* is chosen. The rv’s are then converted to standard 
normal variables u } using equation (59). Lagrange techniques are used to minimize the distance p from 
the limit state curve g=0 to the origin, resulting in the following set of expressions for each of the stan- 
dard normal rv’s: 


Uj =-a-(l, i = l,2,---,n , 

where the gradient unit vector components 0Cj* are defined as 




A 

<9ui 


1 


n f ^ A 2 


S 

j=i 


dg 

V^J J 


(64) 


(65) 


The substitution of equation (64) into equation (59) yields the expression 

— Abij (66) 

Using the assumed values for each primitive rv, these expressions are substituted into the limit state 
function g=0, yielding a single equation for p. An updated design point u * 1 is then obtained by evaluat- 
ing equation (64) using the new value of p. The process is iterated until convergence of p is reached. An 
entire CDF is obtained by using equation (54) to define limit states for the entire range of possible 
response values, and obtaining P for each of these limit states. 

Rackwitz and Feissler continued the development of this method by examining how to develop 
an “equivalent normal rv” for nonnormal rv’s. 28 The starting point of the this procedure is to obtain X* 
using the algorithm discussed above for each rv in X and then determine the actual CDF values, F(x + ), 
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and PDF values, fix*), for that point. An equivalent rv z* is then defined as a function of the design point 
x and the unknown equivalent normal mean and standard deviation parameters, denoted by p x n and c? x n , 
respectively. 


z* = 



(67) 


At this point, F(x*) and f(x ' ) are set equal to a normal CDF and PDF, which are identified by <E> and (p 
respectively, at the same value of x* 


<f>(z *) = F(x *) 


( 68 ) 


<P(z *) 

_n 

°x 


f(x *) 


(69) 


The simultaneous solution of the two equations yields p x n and G x n and is shown graphically in figure 4. 
The new rv with parameters p x n and c? x n is used in the iteration procedure described previously until 
convergence on (3 is obtained. 

Chen and Lind further modified this method by adding a scaling factor y to the matching process, 
calling the overall procedure “fast probability integration” (FPI). 29 This factor is used to multiply the 
normal distribution in order to match not only the CDF and PDF values but also the slope of the PDF 
with that of the actual distribution at the design point, which is usually near the tails. This results in a 
“three-parameter approximation,” where the three equations used to find the equivalent normal distribu- 
tions are: 

y 0(z *) = F(x *) (70) 


7 


~(p{z *) = f(x *) 


(71) 


V7^ 

n~ ( f > ( z *) 


(?f(x *) 

dx 


(72) 


Of course, the new “normal distribution” is not mathematically feasible since the final CDF value will 
equal y. instead of unity, but the shape of the distribution in the vicinity of the design point will be very 
accurate, which is the primary use of this three-parameter distribution. 


Second order reliability methods (SORM), which approximate the limit state with a higher order 
series expansion, have also been examined. The goal of these methods is to describe nonlinear g func- 
tions more accurately to obtain better determinations of the location of the MPP. Rackwitz and Fiessler 
introduced a parameter for measuring reliability in addition to (3 by examining the effect of keeping the 
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second order terms in the Taylor series expansion and considered the different quadratic forms that could 
be used in the analysis. 30 They concluded that FORM was adequate for most engineering problems, 
especially if the rv’s were not too far away from normal. 



The application of FORM and its extensions for nonexplicit limit state functions, as is the case 
for large structural finite element models, requires the use of numerical differentiation to obtain a first 
order approximation of the performance function Y(X), as defined in equation (54). 31 The module FPI 
of the general probabilistic structural analysis code NESSUS implements this differentiation, which is 
defined by the following equation: 

N 

Y(X) = a 0 + £a,(x i -Mxi) , 

i=l (73) 

where 

a o = Y(ju x ) 

N = number of rv’s . (74) 

AY 
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The expansion point a 0 is about the mean value p xi of each of the rv’s. The range of response values for 
a desired probability range is estimated based on the approximate mean and standard deviation as calcu- 
lated using equations (52) and (53). These specific values of y are then substituted in equation (54) to 
generate a limit state equation g(X) for each value. This limit state is used to obtain (3 and the most 
probable point for each desired response value as outlined in the preceding discussion. In addition to the 
errors discussed previously, since the limit state curves are all based upon the expansion at the mean, 
there will be some error for the probability levels obtained if the curvature of the function changes away 
from the mean. 

Alternatively, a CDF can be constructed directly from the approximate mean and standard 
deviation if the distribution is assumed to be Gaussian. This is essentially the original “first order, 
second moment” method as defined by Cornell. For linear functions of normal rv’s, this CDF will be 
identical to that obtained using the MPP iterative search contained in the FORM procedure described 
above. Although it is more direct, implementation of this procedure actually requires an additional step 
since NESSUS/FPI automatically assumes that the input rv’s might be nonnormal and therefore always 
calculates the CDF using the MPP search procedure. Nevertheless, calculation of this CDF outside of the 
NESSUS program using the Gaussian CDF function can prove illustrative in some cases. 

In 1987, Wu and Wirshing developed a procedure for improving the accuracy of probability 
levels obtained by FORM for nonexplicit limit states by using a partial second order expansion. 32 This 
procedure is called the advanced FORM (AFORM). Initially, the FORM is performed and MPP’s ob- 
tained. The limit state is then expanded about each MPP using only the first order and the pure (no 
mixed variable terms) second order terms of equation (47). This is accomplished by performing the 
numerical calculation of g at the mean of the rv’s and at two other points (at least) for every input rv x i? 
resulting in the following new approximation for the limit state: 

n 

g(x) = ^ a i (xj - Xj *) + bj (xj - xj *) 2 ? (75) 

i=l 

where the coefficients aj of the first order terms and bj of the quadratic only terms are obtained by a 
least-squares curve fit of these points. Rather than continue to use this partial second order approxima- 
tion, the quadratic is transformed using algebraic manipulation to a linear function of a set of new rv’s 


n 

g(r) = c 0 + X b i r i , 


i=l 


(76) 


in which 
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4 " b- 

i=l u i 


(77) 
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i )) 


(78) 


The probability transformation rules discussed in Section 3.1 are then applied to x ; to obtain the distribu- 
tions of r i? which will no longer be Gaussian. Therefore, the three parameter normal fit developed by 
Chen and Lind is used to find equivalent normal distributions for q. Using the equivalent normal rv’s, 
the Rackwitz-Feissler algorithm is applied to generate new MPP’s and associated P’s. The new probabil- 
ity value for this P is calculated by multiplying the probability obtained using the normal cumulative 
distribution function of P times the product of the y l scale factor for each transformed rv: 


Pf = 0(y3) P(7i) . (79) 

i=l 

These probabilities have been found to be an improvement over the original values using FORM for 
slightly nonlinear functions. However, for some nonnormal distributions of the input rv’s or for substan- 
tially nonlinear functions, the Rackwitz-Feissler algorithm may not converge to a p and so AFORM is 
not as “robust” as the FORM method. 33 


In cases where the AFORM does not converge or yield realistic answers, use of the “approximate 
statistics” generated at the first stage of the procedure can provide a useful CDF. These statistics are 
generated directly from the partial quadratic form of equations (47) and (48). In addition to the approxi- 
mate mean the approximate median, which is the first term in equation (47) (and which is equal to the 
approximate mean for the first order approximation) is output. These three parameters can be used to 
generate a lognormal distribution X, which is defined as the antilog of a normal distribution Y, with the 
mean p x , median x > and standard deviation o x related to the mean p Y ar| d standard deviation <r Y in the 
following manner: 34 

/n Y = lnX (80) 
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(81) 


This method is accurate for slightly nonlinear functions of normal rv’s, which generally transform to 
lognormal distributions of the new rv’s. 

Wu also developed the advanced mean value (AMV) method, a procedure for updating the 
FORM or AFORM solution by using the original, exact solution for g. 35 Initially, equivalent normal rv’s 
are obtained for any nonnormal rv using the Chen-Lind method. Values of P and X* are then obtained 
for each desired limit state. Each resulting design point is inaccurate, though, for a nonlinear (for 
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FORM) or nonquadratic (for AFORM) limit state, because the approximate limit state function cannot 
track the actual MPP accurately, or if the curvature of the response surface varies away from the mean. 
This will be illustrated in the example in section 3.4. A relatively simple improvement can be obtained 
by transforming the values of the MPP’s to the original nonstandard normal space, and substituting them 
into the exact equation for the limit state, which for structural analysis is usually a finite element model, 
to obtain new response levels associated with the original probability values. Mahadevan explains that 
this incorporates some of the nonlinearity of the actual limit state without having to resort to the com- 
plexities of using a complete second order approximation. 36 This method has been shown to be very 
accurate for several examples given in the literature. 

Further iterations, which are called AMV+, can be performed by expanding the limit state about 
each most probable point instead of about the means as was performed in equation (73). This is a very 
tedious process, and has not been implemented to a significant extent in this research. Both AMV and 
AMV+ have been incorporated in NESSUS, but these automated features were not used in this research 
due to the use of independent FORTRAN codes outside of NESSUS. The incorporation of the automated 
implementation of these methods is one avenue for future work. 

3.4 Comparison of Reliability Methods 

In this section, some of the reliability concepts discussed in section 3.3 are compared and illus- 
trated. Table 1 summarizes the basic premises and areas of applicability of each method. These methods 
are demonstrated with an example. Consider a simple nonlinear limit state function of two rv’s, Xj and 
X 2 , with means p xl and |1 X2 and standard deviations c? xi and g X 2 , respectively. These rv’s are trans- 
formed to standard normal variables x } and x 2 using equation (59). Assume that the exact function of a 
desired response variable Y in terms of the standard normal rv’s, defined as the response surface 
Y(xj,x 2 ), is 


Y(Xj, x 2 ) = 6 + Xj + x 2 + 0.25XJ 2 — 0.25 x 2 2 + 0.25 XjX-, (82) 

as plotted in figure 5(a). The limit state function is therefore formulated according to equation (54) in 
terms of the transformed rv’s, the response function (or surface) Y, and specific response values y: 

g(Xj, x 2 ) = Y(x 1 , x 2 )-y . (83) 

A plot of the function g(x 1? x 2 ) = 0 can be obtained for several different values of y by creating a contour 
plot of the response surface Y, shown in figure 5(b), where the legend identifies the values of y associ- 
ated with the different curves g=0. For a response level of y = 8, the exact CDF probability level can be 
determined using the program “Mathematica.” First, the variable x 1 is expressed as a function of x 2 , 
requiring the solution for the roots of the polynomial g = 0 formulated using equation (83). The root 
corresponding to the g=0 curve at y=8 is 


xj (x 2 ) = -2 - 0.5x 2 -1.11 803^9.6 - 1 .6x 2 + x 2 2 • 


(84) 



Table 1 . Comparison of reliability methods. 


Method 

Method of 
Approximation of 
Limit State 
Function 

Method for Obtaining 
CDF 

Types of rv’s 
Allowed 

Accuracy for Linearity 
of Limit State 

FORM (first order 
reliability method) 

Linear 

approximation 

MPP iterative search 

All types of rv’s 

Robust, medium level of 
accuracy for all levels of 
linearity 

Linear Approx. Statistics 

Linear 

approximation 

Gaussian Distribution 
using p and a 

Normal rv’s only 

Accurate for linear 
functions 

AFORM (advanced 
FORM) 

Linear 

approximation 

MPP iterative search, 
linear transformation of 
rv’s, improved p 

All types of rv’s 

Accurate for slightly 
nonlinear functions 

Partial Quadratic 
Approx. Stats 

Partial quadratic 
approximation 

Lognormal Distribution 
using and x 

Normal rv’s 

Accurate for quadratic 
functions 

AMV Update of FORM 
or AFORM 

Linear or partial 
quadratic 

Update response value 
associated with MPP 
using exact limit state 
calculation 

All types of rv’s 

Improved results for 
nonlinear functions using 
FORM or AFORM 


The double integral of the Gaussian probability density function bounded by the limit state is then 
calculated numerically over the domain -<» < xj < xi(x 2 ) , - °° < x 2 < 


P(Y < 8) 


J Xl ^j 2 ' > exp{^r(-x 1 2 


x 2 2 )} 


2 jt 


dxi 


dx 2 


(85) 


This yields an exact probability of 0.90605. 

The FORM method is applied by initially generating a linear approximation g lin of the function 
Y by expanding about the mean, resulting in 

glin = 6 + x l + x 2 - y • ( 86 ) 

A plot of this surface is shown in figure 6(a). As with the exact function, now termed g exact , plots 
of g lin = 0 for various values of y can be made using a contour plot of the linear response surface 
Y lin = 6 + x j + x 2 , as shown in figure 6(b). To illustrate the MPP searching procedure in FORM and 
some of the errors involved with the linear approximation, the plot of the Gaussian PDF and its contour 
plot are generated (see fig. 7(a)-7(b)). The values of the PDF are identified on the contours. The MPP 
searching algorithm is performed graphically by superimposing all three contour plots in figure 7(c). 

The MPP on g lin = 0 at a specific value of y = 8 is first determined by locating the point on the linear 
contours closest to the origin, specified at point “A.” The value of the CDF is obtained by graphically 
determining the coordinates of the MPP to be (1,1), calculating the distance (3 to the origin, and using the 
Gaussian CDF function to calculate a probability level of 0.921 35 for y=8, an error of +1 .7 percent from 
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(a) Plot of exact response surface Y(x-|, x 2 ). 



X 

(b) Contour plot of response surface showing g=0 curves. 


Figure 5. Plots of exact response function g exact (Xj,x 9 ). 

the exact value calculated above. This example illustrates the error due to the inability of the linear re- 
sponse surface to track the MPP’s accurately, which, for the exact limit state, is approximately at point “B.” 

The AMV method is applied by substituting Xj and x 2 at the MPP “A” into Y exact (Xj,x 2 ) and obtain- 
ing an updated response value of 8.75. This is still associated with the original probability value obained 
using FORM of 0.92135. For comparison with the exact formulation, the AMV procedure can be iterated to 
find the point on the FORM MPP locus (the gradient of the linear response surface) that has an exact y 
value of 8. Since x t and x 2 are equal along this locus, the MPP is solved for in the equation 0=2 x l +0.25 
Xj 2 -2 to obtain point “C” at (0.899,0.899). This results in a probability value of 0.8982, which has an error 
of -0.87 percent from the exact value for y=8, an improvement over the original FORM result. 
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(b) Contour Plot of response surface showing g=0 curves 


Figure 6. Plots of linear approximation of response (g lin x lv x 9 ). 
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(a) Gaussian PDF of two standard normal rv’s. 




gexact(y=7)=0 
9exact(y=8)=0 
| giin(y=7)=o 
j| g lin (y=8)=0 


(c) Superposition of contour plots to obtain MPP’s and probability. 


Figure 7. Obtaining MPPs for FORM, AMY. 
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The AFORM method is illustrated by representing the limit state in a partial quadratic form 
defined as g pquad : 


gpquad = 6 + Xj + x 2 + 0.25xj 2 - 0.25 x 2 2 - y . (87) 

Plots of the partial quadratic response surface and its contour graph are shown in figure 8 (a) and (b), 
respectively. The MPP for y = 8 is estimated graphically to be point A at (1. 163,0.584), and a CDF value 
of 0.9035 is calculated using the Gaussian CDF function, an error of -0.3 percent from the exact value. 
The exact determination of the MPP point for nonlinear limit state functions is performed using the 
Rackwitz-Feissler algorithm described earlier. The contour plot of the exact equation and its estimated 
MPP at point “B” is superimposed for reference. This result shows a considerable improvement using 
AFORM compared with FORM. 

For illustrative purposes, the MPP technique for determining the CDF, which is equivalent to 
integrating the probability volume bounded by a linear slice tangent to the g=0 curve at the MPP, is 
applied for the exact function g exact =0 for y=8. Using the estimated MPP at point “B” of (1.7232, 
0.5514), the distance to the origin is calculated and the Gaussian CDF function yields a probability value 
of 0.88605. This value is less than the exact number obtained previously since the contour g exact =0 for 
y=8 curves away from the mean, yielding a larger probability volume than the linear slice method. Since 
this function is substantially nonlinear, this error is -2.2 percent, which although small, is of the same 
magnitude as the errors due to tracking of the MPP. 

The entire CDF for this function was then examined. The distribution for the exact limit state 
was performed using the “Mathematica” integration methodology described above and verified using a 
40,000 sample MC run in NESSUS, which can accept user-supplied explicit limit state equations as 
input. These results are compared with NESSUS applications of FORM and AFORM in figure 9. Calcu- 
lating the error based on the exact CDF value, the AFORM varies from -0.4 percent error at y=9 to -7.5 
percent error at y=6, while the FORM errors are all slightly larger. 

3.5 Approximate Methods — Perturbation 

The perturbation method is a method for deriving an analytical expression for the change in 
value of a desired response variable as a direct function of the perturbation of the system matrix. For the 
eigenvalue problem, as explained in detail by Meirovitch, 37 it is seen to be a method of obtaining the 
variation in the eigenvalues and eigenvectors of a system due to a change in the mass and stiffness 
matrices without having to resolve the characteristic polynomial. Consider the N dof structural free- 
response problem, which can be expressed as an eigenvalue problem in the form of 

[K]^} 1 =X i [M]{(f)} i , (88) 

where {(f)} 1 is the i’th modal vector. Usually, this can be converted to the standard eigenvalue problem by 
inverting [M] and premultiplying both sides to obtain 

[AKtf-AW . (89) 
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(c) Superposition of contour plots to obtain MPP’s and probability. 


9exact(y=7)=0 
j gexact(y=8)=0 

gpquad(y=7)=0 
| 9pquad(y=8)=0 


Figure 8. Obtaining MPPs for AFORM. 
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Figure 9. Comparison of methods for obtaining CDF of explicit nonlinear limit state. 


where 


[A]=[M]“ 1 [K] . (90) 

Now, let [A] be perturbed by a small variation such that 

[A] p =[A] 0 +[8A] , (91) 


where 


[A] 0 = [M] -1 [K] = original unperturbed system matrix 
[8A] =small change in A (matrix quantity) 

[A p ] = perturbed matrix . 

The eigenvalues and eigenvectors can be expressed in a perturbed form similar to that used above for 
[A]: 


A p= A q+ 5A 1 

(92) 


(93) 


where the subscript “o” represents the original, unperturbed value. 
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The expressions for the perturbed matrices are now substituted into equation (89): 

H + [«])({*}' + W) = (A + SXpi + {*}') . (94) 

This set of N equations in N+l unknowns can be solved for in terms of the known original quantities 
([A 0 ], A 0 , {<j)} 0 ), and the size of the perturbation of the elements in [A] and [8A], by multiplying out 
equation (94) and assuming the effect of the second order terms is negligible. This enables a direct 
expression for the perturbation of A and <|) without resolving the eigenvalue problem. These results are 

a‘= W ‘ T (pA]XC (95) 


■SI®] 1 = 


f wjVM 

jLmi 2 i _ ok 
k=l A o A o 



i,k = l,...,n ; i*k , 


(96) 


where [O] 1 is the complete modal matrix. These first order correction terms can be substituted into 
equations (92) and (93) to obtain the new eigenvalues and eigenvectors. 

In the late 1960’s, this method was applied by Collins and Thomson, 38 Kiefling, 39 and Collins, 
Kennedy and Hart, 40 to probabilistic structural dynamics. Initially, they substituted [K] and [M] back 
into the perturbation equations to obtain the derivative of A with respect to each rv r u : 


dV_ 

dr u 


M 

r( 

d[K] A J <7[M] ^ 

M 

di u 


<2r u 

J 



1 


(97) 


The standard deviation of A! is then obtained by substituting /U for g and using these eigenvalue sensitivi- 
ties in equation (53), which resulted from a first order Taylor series expansion. The mean was obtained 
by directly applying equation (52) for A, equal to g (i.e., calculating the eigenvalue for the mean mass 
and stiffness matrices). The partial derivatives of the elements of the stiffness and mass matrices can be 
obtained by taking the derivative of the analytical expression for each element of the matrix with respect 
to each random property. This method does not take into account nonnormal rv’s, nor does it try to 
describe the eigenvalue distribution at the tails of the cumulative distribution functions since it assumes 
a small perturbation about the mean value. Hart and Hasselman 41 incorporated modal truncation into 
equation (97) so that the entire modal matrix would not have to be obtained a priori. They used compo- 
nent mode synthesis to derive similar analytical expressions for the system eigenvalue matrix [A] as a 
function of the modally reduced substructure stiffness and mass matrices. 


In 1985, Pierre performed an extensive study of the analysis of structural systems with parameter 
uncertainties 42 He formulated the “statistical perturbation method,” which incorporates the correlation 
between input rv’s and response variables (in this case, the eigenvalues), with the perturbation method. 
An expression is derived for the covariance matrix of A as a function of the covariance matrix of {r }, 
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which are the perturbed parameters or input rv’s, using a first order Taylor series approximation: 


[ cov a] = 


' dX~ 

~dA~ 

[cov r ]| 

~dX~ 

"<9Alj 

_<?A_ 

_dr . 

_dA_ 

UJj 


( 98 ) 


For an uncorrelated system, the diagonal elements of [covjJ, which are the standard deviations of the 
response variable, will be identical to the results obtained using equation (53). 

If the chain differentiation for the first two terms in the righthand side of equation (98) is carried 
out, the result will be equal to the eigenvalue sensitivity obtained in equation (97). As with the above 
methods, the mean of X is simply the eigenvalue of the unperturbed system. Pierre also was able to 
improve upon the method described in Meirovitch’s textbook for deriving the partial derivative of X with 
respect to the elements of the system matrix [A] by applying the method developed previously by 
Collins, et al., shown in equation (97): 


dX _ SX' _ 

’ <»> 

where r and s denote the row number and column number of the element in [A], respectively. Obtaining 
the covariance matrix is quite tedious; calculating the partial derivatives of X with respect to each ele- 
ment of [A] is made easier by using the analytical expression of equation (99), but the calculation of the 
partial derivative of each element in [A] with respect to every input rv may have to be performed nu- 
merically, which would be extremely intensive both in terms of central processing unit (CPU) and 
memory requirements. Pierre also derived the covariance matrix using the second order expansion of the 
Taylor series shown in equation (48), but shows that it is “unsuitable for a large number of parameter 
variations.” 

3.6 Comparison of Probabilistic Methods 

Mahadevan has written an excellent unpublished book for NASA/Marshall Space Flight Center 
that reviews and contrasts the various techniques used in the engineering community for evaluating 
probabilistic structures. 43 In this reference, he outlines several approaches for obtaining the limit state 
function defined in section 3.3 for realistic structures modeled using the finite element method. The first 
is the MC method, and it is applied by literally running the finite element model for each sample set of 
input rv’s. This has actually only been made possible in 1997 by an addition to the NESSUS code that 
allows automatic variation of input rv’s in the commercially established finite element code NASTRAN, 
submission of the sample runstream, and statistical processing of the desired response variable. One of 
the achievements of this research was performing a substantial amount of debugging and testing neces- 
sary for using this capability. Even with the new code, a 500-sample run of a realistic bladed-disk struc- 
ture still takes 24 hours to run on a CRAY supercomputer. 

Another method outlined by Mahadevan is “perturbation-sensitivity analysis,” in which an 
approximate closed-form limit state function is obtained. This is done in two ways. The first is the finite 
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difference approach as performed in FORM by equation (73). Mahadevan refers, to this method as the 
“brute-force” method because the response variable is solved directly for a change in every input rv. For 
example, if the response variable is a system eigenvalue for a N-dof system, the steps are as follows: 

1(a) Vary input rv r i by Aq, where i=l,...,d 

1(b) Calculate new matrix [A]=[M] _1 [K], which has N 2 elements (100) 

1 (c) Solve new eigenvalue problem, yielding new X 


< 3/1 

1 (d) Estimate ^ 


^new ^old 
? new — 7 old 


( 101 ) 


1(e) Substitute directly into equation (73). 

The second technique for “perturbation-sensitivity analysis” is the “classical perturbation” 
approach, which uses the chain rule to calculate the partial derivative in equation (101). The steps are as 
follows: 

2(a) and 2(b) Same as 1(a) and 1(b) above 


<3[A] 

2(c) Calculate the partial derivative of each element in [A] with respect to each input rv r, ^ j > 
either numerically or analytically if a closed form expression is available 

2(d) Using equation (99), calculate the partial derivative of X with respect to each element in the 
system matrix [A] 


dX 

2(e) Carry out the chain rule to obtain ^ 


N N 

XX 

r=ls=l 


<9/1 <9A rs 
<9A rs <3rj 


( 102 ) 


2(f) Substitute directly into equation (73). 

Steps 2(d) and 2(e) are equivalent to carrying out equation (97). As mentioned previously, re- 
peated eigensolutions are not necessary for the classical perturbation method; however, the number of 
calculations required is substantially greater than for the direct finite difference approach because of step 
2(c). However, if a closed-form solution for the partial derivatives of the elements [A] with respect to 
the input rv’s is available, this approach would be advantageous. This is why the method has been 
chosen for simplified models of mistuned bladed disks, as will be discussed in chapter 5, where the 
matrix elements have all been expressed explicitly in terms of the input rv’s. For a finite element model, 
these partial derivatives would have to be estimated numerically, so the far fewer number of calculations 
required for FORM would be preferable even though each calculation is more complex. 


37 



This chapter has reviewed the necessary probabilistic groundwork for development of the PDS 
method for analyzing the dynamics of structures. As seen in the next chapter, applying the PDS method 
required a step-by-step implementation of the methodologies presented since no preexisting code could 
accommodate the combination of component mode synthesis with probabilistics. In addition, knowledge 
of the limits of the methodology and ways to extend those limits are essential in evaluating the success 
of the PDS method and identifying sources of improvement. 
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4. DEVELOPMENT OF PROBABILISTIC DYNAMIC SYNTHESIS THEORY 

AND SIMPLE TEST CASE 


4.1 Development of Probabilistic Dynamic Synthesis Methodology 

Enough background material has now been covered to allow a derivation of the PDS methodol- 
ogy. This procedure is a unique combination of the CMS techniques discussed in chapter 2, probabilistic 
reliability methods from section 3.3, and other advanced topics in dynamics and probabilistics, such as 
correlation theory and application of complex loads onto structures. The chapter concludes with a test 
case applying the PDS method to a simple two substructure probabilistic spring-mass system. Much of 
this work appears in reference. 44 

Consider a structure represented by m=l,...,p substructures. Similar to equation (1), the displace- 
ment coordinates for each substructure can be arranged in vectors containing internal and boundary 
dof’s: 



(103) 


Assume that there exists a set of samples of size s for each substructure. Each sample is modally tested 
individually in a configuration such that the interface locations with other substructures are in a free 
condition. For substructure m, sample i, the test will yield eigenvalues { A,} m4 and the eigenvector matrix 
[<J>] m ’i- In addition, the boundary partition of the residual flexibility matrix [Gbb] nv is obtained from the 
measured boundary drive point frequency response functions of the boundary coordinates, as mentioned 
previously. For simplicity, the notation G in this section will refer to only the boundary partition of the 
residual flexibility matrix. These measurements are used to create the stiffness matrix for that substruc- 
ture as shown in equation (34); for free response, only the kept eigenvalues and the boundary partition of 
the kept eigenvectors along with G are required. For ease of analysis, these values are expressed as a 
single vector { x } m,i of dimension Nm defined as 


( h 


m,i 


/t k 

{^bli 

{^b} k 

Gn 


(104) 


l G bb J 
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where A,j is the j’th eigenvalue of substructure m and {<|> b }j is a vector that contains the elements from the 
j’th mass-normalized eigenvector that pertains to the boundary dof’s. If the entire sample i=l,s of sub- 
structure m is tested, { x } m can therefore be defined as a vector composed of elements that are each an rv 
with measured mean and standard deviation. For this research, they have been defined as the set of 
“dynamic” rv’s, in contrast to “primitive” rv’s like material or geometric properties. 

The distribution for these rv’s is assumed to be Gaussian. This assumption is made both because 
it is the most likely distribution for an rv set resulting from manufacturing variations, as is the case for 
structures, and because it is the most convenient for analysis. Techniques for verifying this assumption 
are used in the test case to be discussed below, and a normal distribution did fit the data fairly well. If a 
nonnormal distribution is assumed, the PDS method can still be used, but the rv’s must be transformed 
to equivalent normal distributions using the Chen-Lind procedure described in section 3.3. Because of 
this assumption of normality, equation (59) can now be used to convert { x } m to { x 1 } m , a vector of 
standard normally distributed rv’s (|i=0, o=l) 

There is clearly a large degree of correlation between these dynamic rv’s. This correlation can be 
quantitatively defined by the correlation coefficient p. For two correlated rv’s x and y, 


P = 


Cov(x,y) 

<T X (7y 


(105) 


where the covariance of x and y is 

Cov(x,y) = E[(x- J u x )(y-jUy)] = E(xy)-At x Py . ( 106 ) 

The value of the correlation coefficient ranges from -1 to +1, and can be easily calculated from the 
measured data using the above equations. A value of -1 or +1 means the rv’s are fully inversely corre- 
lated or fully directly correlated, respectively. A value of 0 means the rv’s are uncorrelated. For the 
standard normal variables { x 1 } , it can be seen that the correlation coefficient is equal to the covariance. 
The correlation coefficients of all the dynamic rv’s to each other can be expressed by creating a symmet- 
ric correlation matrix [C], where the first row is the correlation coefficient of the first variable with all of 
the other variables, etc. 

The correlation matrix is used directly to form a set of independent standard normal rv’s { u } m , 
which is required for the application of the reliability methods. The method chosen initially was to 
invoke an orthogonal transformation of {x 1 } 111 using the Cholesky decomposition of [C], which is the 
lower triangular matrix [L] c m , to uncouple the { x 1 } m coordinates, thereby creating {u} m . 45 This can be 
expressed for substructure m as 


{x'} m = [L] c m {u} m . (107) 

This transformation results in a set of rv’s { u } m that have a mean of zero, as with the set { x 1 } m , but a 
standard deviation equal to the eigenvalues of the correlation matrix [C]. Further discussions of the role 
of the correlation matrix will be presented in later chapters. 
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As described in equation (73), a numerical differentiation of the limit state with respect to each 
rv must be performed. This differentiation is performed about the mean value of the limit state. As 
described in section 3.3, an entire CDF of the response variable is obtained by creating a limit state 
function for various response values in the range. If the response variable is chosen to be the first system 
eigenvalue A, the limit state function is therefore 

g(X) = A, sys !(X)-y , 108) 

where the performance function Y(X) is equal to A, sys *(X), formed from the collection of {x} m vectors 

P 

and having length N = ^N m ? anc j y j s a specific eigenvalue in the response probability range. This 

m=l 

eigenvalue problem is solved for the mass and stiffness matrices formulated according to the residual 
flexibility CMS method as shown in equations (33) and (34). 

The numerical differentiation requires the formulation of the substructure system matrices for the 
mean value of each of the elements of { u } m and a formulation for a variation by some small amount of 
each element of { u } m while all the other elements are left at their mean. For the test case discussed in 
the next section, the variation of the rv’s was chosen to be 50 percent of a standard deviation. Since the 
elements of u are standard normal variables, this results in a series of “load” cases, with the first case 
being 



r.5' 


0 

c 

II 

... o 

II 

A 

0 

> 


0 


0 
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where { u } m is of dimension N m . As can be seen, N m vectors { u } are required for m=l,2,...,p substruc- 
tures in the total structure. This set of { u } a , { u } b , ..., { u } p is then back-transformed, first to the corre- 
lated standard normal set { x 1 } m using the transpose of [L] c m , then to the correlated, nonstandard normal, 
original rv’s { x } m using equation (59). The new vectors { A. } , [O b ], and [G] are extracted from the vector 
{ x } for each substructure and used to calculate a stiffness matrix for that substructure using equation 
(34). As shown in equation (33), the mass matrix is generated by simply placing an identity matrix in the 
generalized coordinate partition on the upper left of the complete system matrix. Rather than using zero 
partitions for the rest of the matrix as specified, small dummy values are instead used on the diagonal of 
the boundary partition on the lower right of the mass matrix to provide numerical stability for the 
eigensolution routine. The global system mass and stiffness matrices are now assembled by directly 
coupling the substructure mass and stiffness matrices. This is accomplished by ordering the “kept” dof ’s 
of each substructure sequentially in the system matrices and adding the boundary partitions together. The 
system eigenvalues are then obtained, and the chosen response variable, in this case the first system 
eigenvalue, is isolated. This process is then repeated, with a “load” case for every rv in substructures 1 
through p. For example, the second load case will be identical to the first except that the first element of 
{ u} a will be zero and the second element will be 0.5. After the response variable A sys 1 is obtained for 
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every input rv variation, the numerical partial differentiation of the eigenvalue with respect to each of 
the rv’s in the complete vector of independent rv’s{U}, defined as: 


M = 


M a 

M b 

{u}P 


( 110 ) 


can now be estimated by the first order finite difference technique. 

The actual generation of the limit state function of X sys 1 with respect to each element of {u) is 
performed within the code NESSUS. Each “datapoint,” consisting of a “load” case { U } and its corre- 
sponding value of A.jyg 1 , is copied into the NESSUS/FPI module. In addition, several other parameters 
are established for the NESSUS run. The first parameter specifies whether a linear or quadratic curve fit 
will be used for the limit state function. The next instruction determines whether the FORM, AFORM, 
or another reliability methodology will be used to obtain the most probable points. P’s, and CDF of the 
response variable. Finally, either the set of response values for which a probability level will be obtained 
or the set of desired probability levels for which a corresponding response value is obtained is input. In 
this case, since each element of { u } is only varied by a single value, only a linear approximate limit state 
function can be obtained, so only the FORM can be used. 

The output for this computer run consists of a CDF, p values, and the values of { u } that locate 
the most probable point X* in N-dimensional space for each p value. This solution is also called the 
“mean value” solution because it results from an expansion about the mean of the limit state. The solu- 
tion can now be updated by applying the AMV method explained in section 3.3, which should correct 
some of the inaccuracies due to nonlinearity of the limit state. For this method, each element of the most 
probable point vector {U* } is back-transformed to an original set of rv’s similar to the back transforma- 
tion of the load cases performed earlier. Each of these original sets { x* } is then used in the residual 
flexibility formulation of the substructure stiffness matrices, the substructures are coupled, and a new 
A, sys 1 is calculated. This value is now associated with the previously calculated p value for that most 
probable point. A new CDF is therefore created by shifting the response levels for each of the probability 
levels from the value originally obtained using FORM to the new, updated value. 

4.2 Test Case 

A test case was formulated to further develop and validate the proposed technique. The test case 
consisted of a two substructure, 7-dof spring-mass system (fig. 10). 
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Full Model 



I Sub. a Sub. b 

Figure 10. Test case system. 


To simulate the modal test phase of the procedure, an MC simulation was performed of each substruc- 
ture using 5,000 samples. To achieve complete probabilistic generality, each spring in the system was 
assigned a normal distribution with a p k of 200 and G k of 1 0, and each mass was assigned a normal 
distribution with a ji m of 1.0 and G m of 0.05, yielding coefficients of variation of 0.05 for all the system 
input “primitive” rv’s. As discussed previously, since this initial numerical simulation is performed to 
represent measurements of the dynamic characteristics of a physical population, the distribution of the 
underlying primitive variables is not critical; normal distributions were therefore chosen for conve- 
nience. The MC random vectors were then used to create the mass and stiffness matrices for the sub- 
structures (5,000 for each) and a modal analysis run on each substructure to obtain sample sets for their 
eigenvalues { A.} m and eigenvector matrix [®] m , m=a,b. This analysis was performed using FORTRAN 
and a set of matrix solution subroutines FORMA. 46 

Even if maximum accuracy is desired, only three of the four modes for each substructure can be 
“kept” for the analysis; since the interface locations are modeled in the free condition, this introduces an 
extra mode for each substructure, which if kept, would cause the model to be overspecified. Instead of 
measuring the boundary partition of the 4x4 residual flexibility matrix [G] m using drive point fre- 
quency response measurements, as would be performed for the actual application of this method, it was 
analytically calculated directly from the modes that had been chosen to be truncated as discussed by 
Martinez: 47 


[G] = 



( 111 ) 


where, in this case, k=3 and N=4. It was later discovered that if equation (1 1 1) is chosen as the method 
used to obtain G, substantial error results if all of the deleted modes are not used; this necessitates the 
use of other methods for realistic structures, as will be discussed in chapter 5. 

Proper treatment of the rigid body mode of substructure “b” is also important for the test case. 
Since the first mode is a zero frequency mode by definition, the first eigenvalue is not actually an rv. To 
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handle this in the coding, a mean and standard deviation of zero was manually attached to the variable 
so it could be easily kept as part of the rv set. It is important to note that the mode shape for this mode is 
random; it is a function of the random mass input. By imposing this random condition on the masses, 
complete generality of the problem was maintained. If only stiffness was chosen to be random, as has 
been the case for many other analyses of probabilistic systems, such as with bladed disks, the entire 
effect of the variable rigid body mode shape would be lost. The statistics on the dynamic characteristics 
of each substructure {?i} a , { A, } *\ [‘I*] 3 , [0] b , [G] a , and [G] b and the correlation matrix were then calcu- 
lated. These statistics are listed in table 2. 


Table 2. Statistics of dynamic characteristics. 


Substructure A 

Eigenvalues 

Mean 

Standard Deviation 

1 

30.431 

1.2702 

2 

246.74 

10.336 

3 

552.89 

23.014 

Eigenvectors, boundary location only 

Eigenvector 

Mean 

Standard Deviation 

1 

0.7073 

0.01155 

2 

-0.70725 

0.028546 

3 

0.70688 

0.058827 

Residual Flexibility (one boundary point only) 

Location 

Mean 

Standard Deviation 

1 

6.46E-04 

1.6159E— 04 


Substructure B 

Eigenvalues 

Mean 

Standard Deviation 

1 

0 

0 

2 

150.51 

6.706 

3 

489.03 

21.916 

Eigenvectors, boundary location only 


Eigenvector 

Mean 

Standard Deviation 

1 

0.53466 

6.9105E-03 

2 

0.75592 

2.1945E-02 

3 

-0.55443 

7.4042E-02 

Residual Flexibility (one boundary point only) 

Location 

Mean 

Standard Deviation 

1 

6.46E-04 

1.5347E-04 


The listed quantities comprise the vectors { x } a and { x } b as described in equation (103). A distri- 
bution characterization procedure was also applied to the distributions to verify the normal distribution 
assumption. 48 Partial results of the distribution characterization routine for one of the rv’s are shown in 
table 3. The “W” statistic is a goodness-of-fit test developed by Wirshing and Carlson, where a smaller 
number indicates a closer fit. The results show that the data are well represented by both normal and 
lognormal distributions, with the lognormal being slightly better. The CDF values for each distribution, 
however, indicate that the curves for the two distributions predict almost exactly the same value, which 
can be the case for a particular set of lognormal parameters. The assumption of normality was therefore 
deemed to be accurate. 
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Table 3. Partial results of distribution types routine. 


Substructure B, Eigenvalue 2 

W statistics (all types included) 

Normal: 

0.00948 

Exponential: 

0.30339 

Weibull: 

0.04448 

EVD: 

0.04323 

Lognormal: 

0.00934 

Normal Distribution 

Parameters 

Sample mean = 150.61, 

sample std. dev. = 6.502 

Normal Distribution 

CDF fit to data 

Response Value 

CDF Value 

143 

0.1210 

147 , 

0.2895 

150 

0.4628 

155 

0.7503 

Lognormal Parameters, base e 

mu = 5.014, 

sigma = 0.43201 

Lognormal Distribution CDF fit to data 

Response Value 

CDF Value 

143 

0.1193 

147 

0.2947 

150 

0.4713 

155 

0.7540 


These rv’s were then grouped together for transformation, with the ordering as follows: sub- 
structure “a,” boundary eigenvectors, eigenvalues, and residual flexibility term; then substructure “b” 
boundary eigenvectors, eigenvalues, and residual flexibility term. A matrix composed of cases of {u} 
vectors was generated and multiplied by the correlation Cholesky decomposition matrix [L] c to obtain 
{X 1 }, the set of correlated standard normal rv’s. These were then converted to their nonstandard normal 
distributions and used in the residual flexibility substructure stiffness matrix. The substructures were 
then coupled together and a modal analysis was performed on the system matrices. The first system 
eigenvalue for each case along with its load case, {u}, was then input to the NESSUS/FPI routine. For 
this test case, it was erroneously assumed that the variance of the independent rv’s {u} were equal to one 
instead of the eigenvalues of the correlation matrix; this may have caused some error, although 
Mahadevan has theorized that it may not be significant because of the rapid dropoff of the response 
probability surface away from the mean values of these normalized rv’s. 49 

The FORM solution of X sys 1 from this routine is then used to plot an initial estimate of the CDF. 
The MPP’s for each requested limit state were recorrelated and transformed back, following the AMV 
approach, into the original set {x} m ; this was used to create new substructure stiffness matrices, which 
were coupled to form the global system matrices for eigenanalysis. One of the most probable points is 
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shown in table 4. The updated system fundamental eigenvalues (one for each MPP) were plotted to form 
a new CDF. The comparison of these two CDF’s is shown in figure 1 1 . 


Table 4. Sample MPP output. 


Fundamental eigenvalue response value = 8.101 
P (^)< 8.101 = 0.01 


Most Probable Point 

Substructure A rv 

Value 

Substructure B rv 

Value 

ul 

-0.4863 

ul 

-1.5205 

u2 

-0.1692 

u2 

0.2641 

u3 

-0.05685 

u3 

-0.2361 

u4 

-1.5807 

u4 

0.0 

u5 

-0.1221 

u5 

-0.3057 

u6 

-0.0269 

u6 

-0.0473 

u7 

0.211 

u7 

0.2267 



Figure 1 1 . Comparison of PDS FORM and AMY with CDF’s for seven spring/mass test case. 


There are several choices for verification of the PDS method. For this case, the baseline chosen 
was an MC analysis performed on the same nondeterministic spring-mass system but using the original 
primitive rv’s of mass and stiffness. The system eigenvalue is obtained directly from each 
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unsubstructured sample. This is the full model approach; it yields the closest available approximation to 
the exact solution, but it does not isolate the various approximations used in the PDS method. To isolate 
the effect of substructuring, for instance, it may be more illustrative to perform art MC simulation of the 
two substructures individually and then couple them together to obtain the system eigenvalues. These 
various applications of MC are examined in later sections of this thesis. 

For comparison, the CDF for the full model is superimposed on the previous CDF’s from the 
PDS method in figure 12. A very small amount of error is indicated graphically. To identify the error 
quantitatively, the amount of variation of the fundamental eigenvalue from its mean value at selected 
CDF probability values for the PDS method was compared to the spread for the full model. The result, 
shown in table 5, indicates that the deviations from the mean as computed by the AMV and MC methods 
agree to within 5 percent. In addition, the mean value of the fundamental eigenvalue computed by the 
AMV method is 8.751, which is only 0.5 percent higher than that computed using MC of 8.707, and the 
AMV standard deviation is 0.276, which is only 1.8 percent less than the MC standard deviation of 
0.281. These results were considered acceptable based upon the measure of error used. 



Fundamental Eigenvalue 


Figure 12. Comparison of CDF's using PDS FORM and AMV with MC. 
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Table 5. Dispersion error of PDS versus full MC model. 


CDF 

Value 

AMV 

Delta From 
AMV Mean 

Monte Carlo 

Delta From MC 
Mean 

% Error 
of Deltas 

0.010 

8.116 

-0.635 

8.094 

-0.613 

-3.46 

0.050 

8.300 

-0.451 

8.266 

-0.441 

-2.23 

0.100 

8.399 

-0.352 

8.352 

-0.355 

0.87 

0.159 

8.476 

-0.275 

8.426 

-0.281 

2.22 

0.200 

8.520 

-0.232 

8.467 

-0.240 

3.62 

0.500 

8.751 

0.000 

8.707 

0.000 

Not applicable 

0.800 

8.985 

0.234 

8.944 

0.237 

1.44 

0.841 

9.029 

0.278 

8.987 

0.280 

0.74 

0.900 

9.108 

0.357 

9.072 

0.365 

2.25 

0.950 

9.210 

0.459 

9.169 

0.462 

0.65 

0.990 

9.403 

0.652 

9.344 

0.637 

-2.28 


4.3 Additional Studies of Seven Mass/Spring Problem 

To gain a deeper understanding of the variety of methods available and the approximations and 
errors involved, a variety of additional analyses using this system were then performed using larger cov 
of the input rv’s. It was anticipated that the high cov levels would prove more difficult for the FORM. 
Examination of equation (46) shows that the error in approximating g as a linear function is directly 
related to the size of the second (and higher) order terms in the expansion of g. The size of the second 
order terms is governed by the second partial derivatives of g and by the coefficients of these terms, 
which have the form (x - |i x ) 2 , a measure of the size of the cov. Therefore, any small nonlinearity that 
does exist will be magnified by a large cov. First, an examination of the linearity of the system as a 
function of the independent rv’s {u} was performed numerically. The system fundamental eigenvalue 
was calculated and plottted in figure 13 as an independent function of each independent rv. The range 
for the abscissa is between -2 and +2 times the standard deviation of each rv, which was calculated for 
this examination as the square root of the eigenvalue of the correlation matrix. The legend on the plot of 
the rv’s slul, ..., slu7, s2ul, ..., s2u7 refer to the original order in which the 14 dynamic rv’s were 
placed, with si designating substructure 1 and s2 designating substructure 2. This graph shows that the 
limit state is somewhat nonlinear with respect to several of the rv’s, in particular the first rv for each 
substructure, slul and s2ul, and to s2u4. 

To examine the effect of this nonlinearity, the analysis described in section 4.2 was repeated for 
cov’s of 10 percent for each primitive rv. As expected, the results were found to be slightly worse than 
those obtained for the smaller cov case when using FORM. As with the previous system, the AMV 
update procedure was applied to the most probable points output from the NESSUS/FPI run. The results 
from this did not show any improvement, which may indicate that the AMV method is primarily a 
method of improving the results of nonlinear problems with small cov’s. Further data presented later in 
this section support this hypothesis. Figure 14 compares the results using FORM and the AMV/FORM 
update with the MC baseline analysis for this case. There is some skew at the mean value using the PDS 
reliability methods, but the actual error is only 2 percent. 


48 




Cumulative Distribution Function Fundamental Eigenvalue 



Figure 13. Lineality of fundamental eigenvalue w.r.t. independent rv’s. 



Figure 14. FORM and AMY CDFs versus MC; cov=10 percent. 
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The AFORM method described in section 3.3, which performs a higher order expansion of the 
limit state function, was implemented to try to incorporate the nonlinearities of the system which be- 
come important for the cov = 0.10 case. This required the creation of a minimum of two datapoints of 
variation of each rv in {u} in addition to the value of the mean in order to perform a quadratic curve fit 
Although only three total points are required, in this case four other variations were obtained using the 
following “load” sets, which are placed into a single matrix [U], with the number of rows equaling the 
total number of rv’s in the entire structure N sys = N 1 + N 2 , and the number of columns equaling 
4 N sys + 1 for the mean value set: 

'- 0.2 - 0.1 0.1 0.2 0 ••• 0 " 

0 0 0 0 - 0.2 ••• 0 

[u]= (in: 

0 0 0 0 0 0 


As with the linear methods, this set is recorrelated and back-transformed to a set of original rv’s {x} 
which are used to generate an eigensolution. Each column of [u] and its resulting response variable X. sys 1 
is then input into NESSUS and the “quadratic” g-function option is chosen. This quadratic g-function 
allows the AFORM methodology to be implemented. The results using the quadratic approximation of 
the g-function were noticeably better than for FORM (fig. 15), with an error at the median reduced to 
0.7 percent. The AMV method was also applied to the most probable points output from this quadratic 
run, but the results show almost no improvement over the AFORM results. This may indicate that 
AFORM accounts well for the entire higher-order term in equation (46), including the coefficient, while 
AMY only accounts for excessive nonlinearities expressed by the second partial derivative of g with 
respect to the rv’s. This hypothesis is also supported at the conclusion of this section. 

Since there were approximations and/or possible errors from several different sources for this 
analysis, some attempt at isolation of these sources was made. The first was to run an MC simulation, 
isolating the effect of the reliability methods from the CMS procedure by generating the response vari- 
able using a model generated with the residual flexibility CMS method instead of from a full model, as 
discussed in the last section. As shown in figure 16, which compares the two MC approaches, this 
resulted in a small error from the mean of the full MC of -0.1 percent, indicating that error from the 
above methods is due mostly to the probabilistic approximations involved, not the CMS approximation. 
This was anticipated, since all possible modes were used in the substructure formulation. In addition, an 
MC analysis was run on the full system entirely within NESSUS, as described below, and this solution 
matched very well with that obtained using the FORTRAN codes developed earlier. 

Verification of the reliability methods as implemented in NESSUS was also performed to isolate 
possible sources of error. This was accomplished by directly inputting the full model directly into 
NESSUS using its user-defined FORTRAN subroutine capability, and using the built-in AFORM meth- 
odology to perturb the primitive system parameters (mass and springs), generate the quadratic g-func- 
tion, and obtain a CDF. The results are shown in figure 17, and indicate that the AFORM methodology 
itself performs adequately compared with the MC simulation of the full model. The CDF resulting from 
applying AFORM with the PDS method is shown for reference as well. 
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Cumulative Distribution Function Cumulative Distribution Function 



Fundamental Eigenvalue 

Figure 15. Comparison of AFORM with FORM, AMY; cov=10 percent. 



Figure 16. Comparison of CDF’s for different MC approaches; cov=10 percent. 
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Cumulative Distribution Function Cumulative Distribution Function 
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Figure 17. Verification of NESSUS/FPI AFORM. 



Figure 18. CDF using quadratic “approx, stats;” cov=10 percent. 
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Using the quadratic option also instructs NESSUS to generate the approximate mean, median, 
and standard deviation using the partial quadratic terms from equations (47) and (48) as discussed in 
chapter 3. These parameters are used to create a CDF using an assumed lognormal distribution which 
actually turns out to be closer to the MC results than the AFORM, FORM, or AMV results. The approxi- 
mate statistics result is combined with the NESSUS/FORM and AFORM results in figure 18. 

As seen in figure 15, the AFORM method shows an improvement over the FORM while the 
AMY update of the FORM does not. This difference in performance can be explained by the mechanism 
for inclusion of higher order effects for each method. The AMV partially includes the second order 
information by reassociating values of the CDF with updated response value solutions using the itera- 
tively obtained MPP’s. Since the limit state is expanded about the mean, the MPP for that value is not 
altered but instead stays at {0}, so an update is not possible. The FORM and AMY value at the mean is 
equal to the value of the limit state at the means of the input rv’s (g(X=p x )). Examining equation (47), 
reprinted below for clarity, shows that the error in this value is equal to a term that is directly related to 
the size of the cov times the second partial derivative of g with respect to the input rv, which is only 
nonzero for a nonlinear function: 


/U ♦ /^y ) "b 2 


Vg(fe.A- y ) [ E(x2) _ ^ + ^Yfe.^> [ E ( y2 ) _ ^ 




Even for a function with a small amount of nonlinearity, therefore, the error in the mean can be substan- 
tial if the cov is high, and this error will skew the entire CDF with respect to the MC CDF. Since the 
AMV corrects the value of the response for a given probability, it is therefore most effective in correct- 
ing the shape of the CDF, which is quantified by the standard deviation, rather than a skew, as quantified 
by the mean. The AFORM, on the other hand, uses an approximation of the entire second order term in 
equation (47), so it is able to improve both the orginal response and probability values. Errors in both the 
skew and shape of the CDF can therefore be corrected. 

These studies were also performed on a system with a cov of 0.15. Similar results were obtained, 
and are shown in figure 19. The size of the errors measured at the mean are small, ranging from 
1.6 percent for the PDS AFORM case to about 3 percent for the AMV case, but as seen in the previous 
case, the AFORM is a noticable improvement. 

Both increased computational speed and increased accuracy of input rv’s are among the goals of 
the development of the PDS method. For completeness, comparisons of the speed of the models used in 
the preceding analysis will be presented, but it must be noted that these are not truly direct “apples-to- 
apples” speed comparisons. The bulk of the computational effort in the PDS method was used to convert 
the simulated dynamic data to a set of independent rv’s and then to create the datapoints used to create 
the probability response surface for input to NESSUS. The MC case, on the other hand, was extremely 
useful as an accurate “baseline” case, but directly uses the “primitive” rv’s and so does not have to 
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perform the CPU-intensive step of response surface datapoint generation. In addition, the purpose of this 
test case was to develop the PDS method and to examine the various methods available for obtaining 
accurate statistics of the response variable, not to prove that the method would be faster. Therefore, the 
FORTRAN runstreams developed were not optimized for efficiency and the advantages in the PDS 
method due to the use of the residual flexibility method of CMS were completely nullified by using as 
many dof's in the substructures as are present in the full model to ensure accuracy. 



Nevertheless, the speed comparisons show the PDS method as applied takes the same order of 
magnitude of time to run as the MC procedure. The PDS method uses two main programs — a FOR- 
TRAN code “quad” to create the sample sets of the primitive rv’s, create the load cases, and generate the 
datasets; and a NESSUS/FPI runstream that uses the datasets to create the CDF. The FORTRAN code 
took 80. 144 sec of CPU time and 249 sec of wallclock time on a CRAY T/90 supercomputer. The 
NESSUS run took about 1 sec on an SGI workstation. The 10,000-sample MC simulation of the system 
performed entirely in FORTRAN took 24.8 sec of CPU time and 90.2 seconds of wallclock time. These 
values were in line with expectations based on the reasons stated above. 

In conclusion, for a simple system, even with high values of cov, the PDS method can produce 
good approximations of the CDF using the partial higher order limit state approximation techniques. 
These results were successful enough to verify that the PDS methodology would provide a useful new 
technique for analyzing systems with nondeterministic substructures whose dynamic properties can be 
statistically characterized. It was anticipated, though, that implementation of the methodology on a more 
realistic system would illuminate many more theoretical as well as practical issues that would require 
resolution before PDS could be fully developed. The next chapter, therefore, examines a three-substruc- 
ture system modeled using the finite element technique that brings these issues to light. 
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5. FREE- AND FORCED-RESPONSE ANALYSIS OF REALISTIC SYSTEM 


This chapter describes the analysis of a realistic three-substructure probabilistic system using the 
PDS method. The purpose of this analysis is to identify and resolve problems associated with implemen- 
tation of the method on a realistic engineering structure. Section 5.1 defines the problem motivating the 
dynamic analysis for a nondeterministic structure undergoing dynamic excitation. This definition is used 
to identify a metric for evaluation of the probabilistic methods applied in section 5.2. The deterministic 
and stochastic geometric and material properties of the structure are then described in section 5.3. The 
numerical simulation of the modal testing of the probabilistic substructures, which is the first step of the 
PDS process, is explained in section 5.4. Section 5.5 applies the methodology to obtain the partial free- 
response solution of the structure, which is a CDF of the third eigenvalue of the system. Finally, the 
forced response solution in the frequency domain is presented in sections 5.6 and 5.7. 

5.1 Problem Definition 

It is important to define the problem that is faced by designers and analysts of structures undergo- 
ing dynamic excitation. Frequently, the approach used in industry is to design structures initially to avoid 
resonance with potential harmonic excitation mechanisms, and then if that resonance cannot be avoided, 
to perform a forced-response analysis to evaluate the level of response. The problem statement therefore 
is divided into a free-response and a forced-response stage. For the free -response stage, standard practice 
has been to impose a rather arbitrary 5- or 10-percent band around the natural frequencies. This may have 
some basis in modal testing of the structure (e.g., modal testing of individual turbine blades), but it is not 
accurate if there is coupling between substructures nor is it tied quantitatively to a probability value. 

The probabilistic approach to the free-response problem is to determine the actual “3-c” range of 
the structural natural frequencies about the mean and design the excitation to avoid that range. This 
entails calculating the CDF for each mode and using the natural frequencies at 0.14 percent and 99.86 
percent (a ±3-a range for a normal distribution) as bounds. Therefore, a probabilistic method should be 
evaluated by its ability to yield the correct natural frequency for these given CDF values. In this research, 
1 percent and 99 percent values were the closest probabilities to these 3-c? values that were obtained and 
so will be used instead. 

If this frequency range cannot be avoided by the excitation, which is a frequent condition, then the 
forced response problem must be solved. In this case, standard practice has been to perform a determinis- 
tic analysis using the nominal geometry of the structure to obtain a maximum displacement or stress, and 
then to impose a “safety factor” on the result. This factor is purely experience-based since there is no 
current methodology in use by industry to evaluate the statistical variation of the response. The probabi- 
listic approach for the problem is to determine the stress CDF of the maximum responding location and 
to design the structure to withstand the positive 3-0 response value. For the forced-response problem, 
therefore, a probabilistic method should be evaluated by its ability to obtain an accurate response value 
for a CDF value of 99.86 percent. 
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5.2 System Description 


In order to examine the applicability of the PDS method for a realistic design problem, a struc- 
tural system had to be chosen that met several criteria. First, the structure had to be modeled using 
standard methodology. This was satisfied by using the most widely used commercial finite element 
code, NASTRAN. Next, since tying CMS together with probabilistic methods is an important goal of 
this research, the size of the model had to be large enough to be able to realize a substantial reduction in 
computer time due to dynamic reduction. This requirement is satisfied by using a structural system 
composed of a “disk,” which is made up of 630 quad4 plate elements and constrained at the center, and 
two “blades,” which are each composed of 24 quad4 plate elements (see fig. 20). For ease of analysis, 
the standard assumption that in-plane translation and rotations are small has been applied by constrain- 
ing out these dof’s, thereby halving the size of the model from 4,142 dof ’s to 2,076 dof’s. To pave the 
way for the detailed bladed-disk analysis presented in chapter 6, the disk was assumed to be determinis- 
tic and the blades nondeterministic. To make the structure somewhat generic, a different mean thickness 
was defined for each of the blades. An additional goal requiring a model of this size is that it necessitates 
the use of model reduction techniques for coupling large, realistic substructures that are not needed for 
very small models like the one used in chapter 4. 

Finally, the structure had to possess a “generic” level of randomness not easily defined by a 
single variation in a material or geometric property. A single rv would not be able to capture variations 
in mode shapes that are independent of variations of natural frequency, for instance. To achieve this 
goal, each blade was separated into three sections, with two of the sections having a thickness set to be 
an independent primitive rv. In addition, the density of each blade across all three sections was defined 
as an independent rv, thus giving each blade three independent rv’s. The different section properties for 
the two-bladed disk are shown in figure 20, and details on the geometry and other properties of the 
system model are shown in table 6. This level of randomness is similar to imposing a random field, 50 
which is a methodology for imposing a correlated random spatial field along a structure. The use of this 
discipline is beyond the scope of this thesis. 
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Figure 20. Two-Bladed Disk 
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Table 6. Three-substructure system (“two-bladed disk”) information. 


Disk Substructure — Deterministic 

Mumber elements: 630 
Number nodes: 631 
Number dot’s: 1,893 
Number unconstrained dofs: 1,884 

Boundary conditions: Fixed at center and two adjacent nodes off center, node 2 and 12 
Material: Steel, Young’s Modulus E= 3.0*10 7 , Poisson’s Ratio v=0.29, 

Density p=0.000769 lb-sec 2 /in. 

Geometry: Diameter = 2.51 in. 

Thickness t = 0.1 in. 

Blade A 

Number elements: 24 

Number nodes: 36 

Number dot’s: 108 

Boundary conditions: Free-Free 

Material: Steel, E= 3.0* 10 7 , v=0.29 

p - N^.eS’lO -4 lb-sec 2 /in., 7.69*10 -5 lb-sec 2 /in.) — independent rv 

Geometry: Length = 2 in. 

Thickness: Section 2: t = 0.1 in. (deterministic) 

Section 3: t - N(0.1 in., 0.01 in.) — independent rv 
Section 4: t - N(0.1 in., 0.01 in.) — independent rv 

Blade B 

Number elements: 24 

Number nodes: 36 

Number dofs: 108 

Boundary conditions: Free-Free 

Material : Steel, E= 3.0*1 0 7 , v=0.29 

p ~ N(7.69*10‘ 4 lb-sec 2 /in., 7.69*1 O' 5 lb-sec 2 /in.) — independent rv 

Geometry: Length = 2 inches 

Thickness: Section 5: t = 0.2 in. (deterministic) 

Section 6: t ~ N(0.2 in., 0.02 in.) — independent rv 
Section 7: t ~ N(0.2 in., 0.02 in.) — independent rv 

Unsubstructured System 

Number elements: 678 
Number nodes: 695 
Number dofs: 2,085 
Number unconstrained dofs: 2,076 

Boundary conditions: Fixed at center and two adjacent nodes off center, node 2 and 12 
Other properties same as above; all rv’s independent 


5.3 Monte Carlo Analysis 

To provide a measure against which to measure the accuracy of the methodologies developed in 
this thesis, a baseline had to be generated. As in chapter 4, an MC analysis of the original system with its 
primitive rv’s was chosen to be this baseline. In addition, MC analysis was required in this research to 
simulate the modal testing phase of the PDS methodology to obtain the statistics of the dynamic rv’s. 
However, probabilistic analysis of structures with this level of detail and randomness has not been 
reported in the literature, even using MC, and is one of the unique aspects of this research. One of the 
reasons that so little research and analysis has been done with this type of realistic structure is that the 
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tools have essentially been unavailable. The finite element code contained within NESSUS was capable 
of modeling a structure similar to that described; however, the accuracy of the elements was suspect, and 
the code was developed with very basic research in mind and so was lacking many capabilities required 
for PDS, such as mass normalization of modes. In 1995—1996, NASA/Marshall Space Flight Center 
contracted Southwest Research Institute to develop an interface with NASTRAN. 51 This interface 
allows any geometric or material property to be specified as an rv with a given distribution type and 
associated descriptive parameters (mean, standard deviation, etc.). The finite element model is read 
directly into NESSUS where the probabilistic algorithms described in chapter 3, including MC, FORM, 
and AFORM, can be performed. In addition, an extensive user-defined interface is available, allowing 
implementation of user-created FORTRAN or C subroutines. The PDS analysis of the three-substructure 
model described above was essentially a “beta” test for the NESSUS/NASTRAN interface, requiring the 
resolution of a number of bugs and software usability issues. The interface has now been made a stan- 
dard part of the NESSUS code, which is commercially available. In addition to the unavailability of a 
code for performing the analysis, the supercomputing capability required to run several hundred (at a 
minimum) sequential finite element runs has not been available until the last several years. The runs 
performed for this analysis took up to 27 hr on a CRAY T-90, one of the most modem supercomputers 
available. 

The new NESSUS/NASTRAN interface manages the simulation, both by creating the random 
vector set for the given input rv’s and by automatically submitting the jobs. The actual implementation 
of this automatic submission, however, required the solution of many network and software problems, as 
this was the first real test case for the code. Reliable solutions were eventually obtained, however, for 
use in the research. 

The MC simulation is extremely time consuming, so a much smaller number of runs was used 
for the two-bladed disk system than for the test case discussed in chapter 4. Applying equation (45), a 
confidence level of 98.2 percent is obtained for the 0.01 probability level using 400 samples, which 
therefore should be sufficiently accurate. For verification, 1 ,000 sample cases were also run, and the 
resulting CDF’s matched the 400 sample cases almost identically. The results using 1 ,000 samples were 
used for the analysis of the two-bladed disk, but based upon this agreement, only 400 samples were used 
for the larger bladed disk analyzed in chapter 6. 

5.4 PDS Methodology — Simulation of Statistics from Modal Test 

As with the simple system described in chapter 4, the first step in the analysis is to perform MC 
simulations of the nondetermini Stic substructures to produce the dynamic rv sample set. This sample set, 
if actually measured from a population of substructures, would be significantly more accurate in repre- 
senting the random substructures than using a set of assumed primitive rv’s. One of the first new issues 
that arose for this realistic system was producing a statistically consistent set of rigid body modes for the 
two free-free blades. Since there are 3 dof’s per node that are free, each blade will have three unique, 
zero-frequency rigid body modes. In general though, the shape of the modes (i.e., the combination of 
rigid-body translation and rotation) and the order that these modes are produced by a numerical proce- 
dure is arbitrary, since they are all at the same frequency. The PDS method requires that statistics be 
taken on the mass-normalized mode shape, so a method of producing consistent, reproducible output 
was necessary. 
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Because of this requirement, it was seen early in the analysis that the finite element capability 
built into NESSUS itself was unacceptable since it would not be able to provide this information di- 
rectly. Fortunately, detailed study into NASTRAN capabilities resulted in the discovery that the use of 
the SUPORT card with the modified Householder method of eigenvalue extraction would yield exactly 
what was needed. This solution would produce rigid body modes parallel to the unconstrained coordi- 
nate frame, and the order of these modes would be consistent for the random perturbations in the model 
resulting from the probabilistic analysis. Sample sets for these mode shapes could therefore be obtained 
where, for example, the first mode would always be a translation parallel to the global z axis, the second 
mode would be a pure rotation about the x-axis, and the third mode would be a pure rotation about the y- 
axis. This enabled the generation of dynamic rv’s from these modes with consistent means and standard 
deviations that could be used in PDS. 

Another aspect of numerical eigenvalue extraction methods is that the signs of the modes are 
arbitrary and can be reversed for slight variations of the model. Since gathering the statistics of the 
modes require that they all be of one sign, a subroutine was written to compare the sign of the modes for 
the second and later samples with that of the first, and, if it was different, to multiply all the values in the 
mode by -1. 

Calculation of the residual flexibility matrix from the modal data, which was performed using 
FORTRAN subroutines linked to the NESSUS run, was also different than for the simple case. More 
detailed research and some test runs indicated that significant error would result if all of the truncated 
modes were not used in the calculation using equation (111). Since obtaining all the modes is generally 
not realistic, alternate methods for obtaining the residual flexibility had to be obtained. For the con- 
strained deterministic disk, the system stiffness matrix is first inverted to create the total system flexibil- 
ity matrix, as expressed in equation (16), and the residual flexibility is then obtained by subtracting the 
flexibility of the retained modes, which are, in general, much smaller in number than the truncated 
modes, from the total flexibility according to equation (22). 

For the free-free probabilistic blades, the system is unconstrained so the stiffness matrix will be 
singular. The “inertia-relief’ method, initially developed by Craig, 52 is therefore used to obtain the 
residual flexibility. First define the transformation matrix P according to: 

P = I - MO r <E>J , (113) 


where O r contains the rigid body modes and M is the substructure mass matrix. Now choose any stati- 
cally determinate subset (in this case, of size three), denoted E, from the boundary dof ’s, and partition 
the stiffness matrix accordingly into the E set and the other set, denoted by a subscript O: 


[K] - 


Kee 

Kqe 


Keo 
Koo. • 


(114) 
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Next, consider a flexibility matrix composed of only the inverse of the K EE partition and defined as G s : 


G 


S 


O' 

0 o • 


(115) 


In terms of P and G s , Craig shows that the flexibility matrix for the elastic modes is given by 


q elastic P^G^P 


(116) 


Equation (22) is now applied to obtain the residual flexibility: 


residual _ r ' elastic ^ retained _ ^ e 


0CIclblIC _ 0 


= G'-[4> N ][A N ]- 1 [4> N ]' , 


(117) 


where <I> N are the retained elastic modes. The calculation of these residual flexibility values were calcu- 
lated for each sample during the MC run, and saved so that the statistics could later be generated. 

During this formulation of the analysis, it became evident that the size of the dynamic rv set 
would become intractable for a realistic problem. For this problem, which has only two probabilistic 
substructures and where 20 modes were chosen to be retained (out of a possible 108) per substructure, 
the number of dynamic rv’s would be obtained as follows: 

Number dynamic rv’s = 2 substructures*! 17 elastic eigenvalues 
+ (4 boundary nodes)*(3 dof ’s/node)*(20 retained modes) 

+ [(4 boundary nodes)*(3 residual flexibility dof ’s/node)] 2 } = 802 

It is evident that this size is highly dependent on the number of retained modes. Standard practice is to 
include up to the mode that has twice the natural frequency as the frequency of interest. In this case, the 
twentieth mode for the disk is at 7,596 Hz, and the blades are both above 50 khz. The frequency of 
interest, which was chosen for the third system mode, is only 456 hz, so the chosen modal representation 
should be sufficient. This rule of thumb can be inaccurate in certain cases. Exploring its applicability to 
this case is an avenue of future research. 


Some attempt was made to incorporate this large number of rv’s in the problem, which requires 
alteration of the NESSUS code as well as massive matrix storage issues. Rather than continue along this 
course of action, several assumptions were made to drastically reduce the number of dynamic rv’s. The 
first was to assume that the limit state function was insensitive to the variation in the rotational dof ’s in 
the modes. This reduced the size of eigenvector rv’s from 240 to 80 per substructure. The second was to 
assume the limit state was insensitive to not only the rotational dof ’s but also to the off-diagonal terms in 
the boundary residual flexibility matrix, which reduces the size of that contribution from 144 to 4 per 
substructure. These assumptions do not remove these variables from the formulation of the substructure 
stiffness matrices; instead, it allows the use of the mean value of those variables (or median, as will be 
discussed later). Admire and Tinker 53 examined the effect of completely removing the off-diagonal G 
elements, and calculated a natural frequency error of less than <5 percent between the exact value and 
the value obtained using the RF formulation. 
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After the entire sample set for the dynamic rv’s was created, the program is then directed to a 
subroutine to calculate the correlation matrix for the blades. Recall from section 4.1 that the Cholesky 
decomposition technique is used to decorrelate the rv’s. For this problem, the Cholesky decomposition 
broke down because the correlation matrix was positive semidefinite; for reasons discussed below, many 
of the eigenvalues were almost equal to zero. Due to the ill-conditioning of the correlation matrix, an 
alternate procedure, which used the eigenvalues { A, c } and eigenvectors [<f> c ] of the correlation matrix 
was implemented. The matrix of eigenvectors [<f> c ] was used to construct an orthogonal transformation 
from the independent rv “load set” {u} to the standard normal variables { x* } , which are obtained from 
the original vector of dynamic rv’s { x } using the mean and standard deviations of the elements, as 
performed for the benchmark case. This transformation is shown as follows: 

{ x 1 } = [<3> c ] { u } . (118) 

The set { u } consists of independent rv’s whose variances (square of the standard deviation) are elements 
in {A, c }. 

To verify that the calculation of the correlation was done correctly, in this case the back transfor- 
mation of equation (118) was performed immediately to create {u}and the correlation between these rv’s 
was then calculated. If the elements of {u} were uncorrelated, the correlation matrix would be an iden- 
tity matrix. The result was somewhat different, but was very illuminating. For the rows/columns of the 
new correlation matrix [C u ] that corresponded to eigenvalues on the order of one, the diagonal term was 
equal to 1 .0 and the other elements in the associated row/column were very small. However, many of the 
eigenvalues turned out to be very small and the analogous rows and columns in [C u ] were not diagonal; 
the diagonal values were 1 .0 but the other terms in the row/column were all on the order of one (between 
zero and one) instead of being close to zero. This can be explained by realizing that the actual number of 
independent rv’s in the system is equal to the original number of independent primitive rv’s, not the 
number of dynamic rv’s. Therefore, the rv’s associated with the small X c values are actually insignificant 
and the magnitudes of these values are mainly due to numerical roundoff. The high off-diagonal values 
in the new correlation matrix in these rows/columns therefore indicates a large amount of correlation 
between these numerical artifact “rv’s.” 

The above result can be used to further decrease the size of the rv set { u } . A sum of the eigenval- 
ues of the correlation matrix was calculated, and any eigenvalue < 3 percent of the sum was deemed 
insignificant and the associated rv in { u } was neglected. The value of 3 percent was reached by starting 
at 5 percent and decreasing the cutoff value (thereby including more rv’s) until the final result stabilized. 
This reduced the eigenvalue set from 101 to 1 1. It is unclear why this value did not actually reduce to 
three, which is the number of primitive rv’s per blade, but nevertheless, this reduction greatly facilitated 
the analysis. 
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5.5 Use of Dynamic Random Variables for PDS 


As described in chapter 1, one the main goals of this research is to use the statistics of the dy- 
namic characteristics of the substructures as the rv’s in a probabilistic analysis rather than primitive rv’s. 
These dynamic rv’s can be applied either using an MC approach or using the reliability method ap- 
proach. The MC should provide more accurate answers, while the reliability method will be more 
computationally efficient. 

5.5.1 MC Using Dynamic Random Variables (PDS MC) 

Using the assumption that the dynamic rv’s followed a normal distribution, a vector of each rv 
was created using the data obtained in section 5.4 to generate the matrix [u], where the number of rows 
equals the number of reduced independent rv’s (2*ll/blade=22) and the number of columns is equal to 
the number of samples. This matrix is transformed to the correlated standard normal variable matrix [x 1 ] 
(of dimension equal to number of dynamic rv’s by number of samples) by using a reduced matrix [3> c ] 
(of dimension equal to number of dynamic rv’s in [x 1 ] by number of reduced independent rv’s). Equation 
(118) therefore becomes 


[x 1 ] = [O c ][u] . (119) 

As before, the statistics from the simulated modal testing are then used to transform [x 1 ] to [x], the set of 
original dynamic rv’s for each substructure. Instead of using the mean values for the transformation from 
standard normal rv’s to the original rv’s, the median was used. This is explained in the next section. The 
matrix [x] is then used to create the substructure stiffness matrices which are coupled together with the 
deterministic disk substructure to create a system stiffness matrix for each sample for eigenanalysis. The 
procedure used was essentially the same as for the benchmark case except that the complexity was 
considerably greater. This increased complexity was first evident in the ordering procedure used to 
rearrange the substructures into internal and boundary dof’s and in the method used to tie the boundary 
dof’s of adjacent substructures together. These procedures required the use of integer vector “maps,” 
which label each substructure dof with its total structural dof. This methodology is used widely through- 
out industry for coupled loads analysis. 

A modal analysis is then performed on each of these sample systems to generate CDF’s for any 
particular free-response characteristic. For the purpose of this study, the third natural frequency of the 
combined system was selected. The comparison with the MC baseline is shown in figure 21. The results 
are mixed; from the design problem approach, the errors as seen in the tabular portion of figure 21 for 
the .01 and .99 CDF levels are < 2 percent, so the values could be used with confidence. From a theoreti- 
cal standpoint, the curves do not line up extremely well; this is indicated by the error in the standard 
deviation of about 20 percent. Based on the results of chapter 4, it was anticipated that the performance 
of this method would be better. The potential sources of error are as follows: 
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(1) Transformation error of the statistics of the original primitive rv’s due to inaccuracy of the 

MC simulation of the dynamic rv’s (possibly insufficient number of samples). 

(2) Transformation error in the statistics and recorrelation of the rv’s due to their assumed nor 

normality. 

(3) Error due to truncation of the number of independent rv’s. 

(4) Truncation error due to using the RF formulation of CMS. 

(5) Programming error. 



5.5.2 Using Dynamic Random Variables in Reliability Approach (PDS Reliability) 

The data obtained in section 5.4 are now used to generate the different “load” cases of [u] for use 
in the reliability method of probabilistic analysis as opposed to for an MC simulation. Matrix [u] is 
created in this approach by varying each of its elements by ± 0.3 of a standard deviation, which for that 
rv is the square root of the correlation matrix eigenvalue, as described in section 5.4. The number of 
rows in [u] is equal to the number of probabilistic substructures times the reduced number of indepen- 
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dent rv’s, and the number of columns is equal to this number times two, plus one additional column for 
the mean value set. The superscripts j =1,2 on the matrix variables below refer to the number of the 
probabilistic substructure, and the subscript i =1,...,11 refers to the number of the independent rv: 


-0.3* cr i 

u i 


[u] = 


0 

0 


0.3 *cr 

0 

0 


0 0 0 0 

0 0 0 

0 -0.3 * o ■ 0.3*0' . 0 

u? „J 


( 120 ) 


The response value of interest, the third eigenvalue in this case, was then appended to each column of [u] 
from equation (120) to create a NESSUS/FPI “datapoint.” 

With these datapoints, the AFORM quadratic solution method can be applied by NESSUS. By 
simply changing the first variation in the above approximation to a value very close to zero times the 
standard deviation, the three points can also be used to create a FORM linear approximation. The qua- 
dratic and linear methods yielded virtually identical CDF’s. This implies that the system is very linear 
and that the coefficient of variation of 10 percent is not large enough to magnify the small nonlinearities. 
The linear (FORM) CDF is plotted along with the MC baseline solution of the total system in figure 22. 
Because of the large difference in the CDF’s as seen visually and as indicated by the size of the standard 
deviation, sources of error not previously examined were investigated. 

One subtle but important error found was in the application of the transformation to standard 
normal variables. For the MC baseline solution, the results are obtained by creating a normal distribution 
about the mean values of the primitive rv’s; this distribution is symmetric, so the median is equal to the 
mean. Since the finite element solution is slightly nonlinear as a function of the input rv’s, the mean of 
the solution will not equal the median (see equation (47)), but the median solution will result from using 
the mean (equal to median) primitive rv’s. 

As PDS is implemented for these analytical simulations, there is an intermediate step that intro- 
duces error. Symmetric Gaussian MC distributions are created about the means of the dynamic rv’s, 
which will not be the medians since a slightly nonlinear eigensolution has been performed to obtain 
them. This results in a skew of the entire CDF curve since the method assumes that the dynamic rv’s do 
in fact follow a symmetric distribution. The median of the results (indicated by the response value at a 
CDF level of 0.5) should result from the median of the input rv’s. However, the deterministic value for 
mode three obtained for the medians of the primitive input rv’s using the RF CMS method is equal to 
455.3 Hz, which is 1.4 Hz less than the median value obtained using AFORM of 456.7 Hz. 

It should be noted at this point that the RF CMS method was used for this deterministic calcula- 
tion to provide a consistent basis for comparison with the reliability methods, which also use the RF 
CMS method. This calculation can also serve as a checkout of the RF CMS method itself by comparing it 
with the deterministic value obtained using the original NASTRAN unreduced model using the medians 
of the input primitive rv’s. The NASTRAN value of 456.035 is only 0.735 hz, or 0.16 percent, greater 
than the RF CMS deterministic value. These values are also shown in figure 23. 
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Figure 22. Two-bladed disk, mode 3 natural frequency comparison of MC baseline with 
PDS FORM with mean. 


1 ■ 

0.9' 

«= 0.8 

0 

| 0.7 

.1 06 

| 0.5 

° 0.4 

1 0.3 
E 

o 0.2 
0.1 
0 

438 442 446 450 454 458 462 466 470 474 

Third Natural Frequency (Hz) 




1 II 1 1 1 1 1 1 1 1 1 1 























FORM Using Means 
NASTRAN Median Determ. 





A 






£ 

" 


J 

4 















k RF CMS Median Determ. 

— * — FORM Using Medians 





/ 























It 

/ 



















MC Baseline (1 ,000 sample) 



f 














































/ 


I 

































/ 

J 


/ 

\ 































y- 

/ 

/ 

. , 

v 






























/ 

y' 





























































Parameter 

MC Baseline 
1.000 

FORM Using 
Mean 

% Error 

FORM Using 
Median 

% Error 

Mean 

456.1482 

456.725 

0.13 

455.297 

-0.19 

Standard Dev. 

5.497802 

2.25244 

-59.03 

7.08748 

28.91 

Median 

456 1161 

456.725 

0.13 

455.297 

-0.18 


Figure 23. Two-bladed disk mode 3 natural frequency comparing FORM using median 
versus mean. 


The error resulting from using the means can be corrected by simply ensuring that the median 
primitive values are carried through the procedure to generate the median results. This is accomplished 
by using the medians xof the calculated dynamic rv’s as the “mean” in the transformation to standard 
normal coordinates: 


{x} = |x l |(7+x . (121) 

Although there is still some error inherent in the normal assumption, this correction allows for a direct 
“apples-to-apples” comparison with MC baseline simulations. The results for this correction are shown 
in figure 23 and are clearly better. 

One unexpected benefit from using this correction was the vast improvement not only in the 
skew of the CDF curves at the median, but also in the standard deviation. This change can be understood 
by examining the variation of the response variable to the input rv’s when the expansion is performed 
about the medians instead of the means. The expansion is shifted to an entirely new location in the 
hyper-dimensional response surface space, and the curvature is different. This is indicated by looking at 
the variation of the response variable with respect to a single independent rv, ul6, for an expansion 
about the median and about the median. As seen in figure 24, the sensitivity curve is much steeper when 
the expansion is about the median than the mean. This indicates a higher standard deviation, or spread in 
the CDF curve, since the same variation in the input rv will yield a larger change in the response vari- 
able. This increased standard deviation is seen in the final comparison of results in figure 25, which 
shows a much better fit with the MC. The results are discussed in more detail in section 5.5.3. The 
higher order reliability methods were also applied to this problem using the medians instead of the 
means of the dynamic rv’s. Neither the AMV or the AFORM showed any improvement over the linear 
case; this may be due to the very small nonlinearity of the limit state. 

In addition to the methods discussed above, the approximate statistics generated as a first step of 
the AFORM solution can be used to create a CDF, as mentioned in chapter 4. Using the mean and 
standard deviation from the FORM solution with an assumed Gaussian distribution yields a CDF identi- 
cal to the FORM solution, as expected. Since the AFORM solution approximate statistics yield a differ- 
ent value for the median and mean (since a partial quadratic calculation of the mean is used), it is more 
appropriate to use these in an assumed lognormal distribution. The resulting CDF was also virtually 
identical to the linear solution, indicating a strong degree of linearity in the system. 

5.5.3 Results Using PDS Methods 

The results using the PDS MC, PDS reliability, and PDS “approximate statistics” methods are 
compared with the MC baseline CDF in figure 25. The graph can be used to visually compare the total 
match of the CDF’s, while the table should be used to evaluate the errors for the locations of interest in 
the design problem. As with the PDS MC, the actual error for the design values is very small when using 
the PDS reliability methods, so they could be used with confidence, but the discrepancy in the shapes of 
the CDF’s, as indicated by the errors in the standard deviations as well as visually, suggest that further 
improvements could be made. The potential sources of error for the reliability methods are the same as 
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ul 6 percent of Standard Deviation 


Figure 24. Two-bladed disk, mode 3 frequency as a function of ul6, using medians and 
mean values of dynamic rv’s. 



Parameter 

MC Baseline 
1,000 

FORM Using 
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% Error 
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Carlo 

% Error 

Lognormal 

Dist. 

% Error 

Mean 

456.1482 

455.297 

-0.19 

456.601 

0.10 

458.747 

0.57 

Standard Dev. 

5.497802 

7.08748 

28.91 

6.8695 

19.97 

7.231 

23.97 

Median 

456.1151 

455.297 

-0.18 

456 

-0.03 

455.296 

-0.18 

1% CDF Value 

444.6 

438.81 

-1.30 

450.75 

1.36 

438.95 

-1.29 

99% CDF Value 

468.57 

471.79 

Q£2_ 

460.35 

-1.79 

472.35 .. 

0-80 


Figure 25. Two-bladed disk, mode 3 natural frequency comparing MC with PDS. 
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for the PDS MC method with the addition of other inaccuracies as discussed in chapter 3. 

It is interesting to compare computational times at this point. All the analyses were performed on a 
CRAY T-90 supercomputer. The MC baseline analysis took 2,527 CPU sec and 20.3 hr of wallclock time to 
run. PDS MC took 389.5 sec of CPU time and 1,1 12.3 sec (18.5 min) wallclock time, a drastic reduction 
enabled by the continuous generation, storage, and analysis of modally reduced system matrices rather than 
the complete model generation and analysis necessary for each of the MC baseline samples. The FORM 
method took only 20.1 sec of CPU and 64.2 sec of wallclock to create the datapoints for input to NESSUS/ 
FPI, and about 30 sec wallclock to run FPI to perform the FORM solution. This further reduction is due to 
the substantial decrease in the number of solutions generated. As anticipated, both the use of CMS and 
probabilistic methods drastically decreased the amount of time necessary to run the analysis, which, espe- 
cially for a more detailed model, could be a prerequisite for use in design. 

5.6 Forced Response — Deterministic Case 

The use of PDS was now expanded to forced response. In particular, a frequency response solution 
of the system was derived since this type of excitation would prove most applicable to the bladed-disk 
problem. This requires a rederivation of the basic equations derived in chapter 2 for the RF CMS method. 
Initially, the dof ’s of a substructure are partitioned into three sections: x 0 , which are internal dof ’s with no 
external load applied; Xj, internal dof ’s with an external load applied; and x b , boundary dof’s which may or 
may not have external load applied. The equation of motion is therefore 


[M]< 

V 

x i 

> + [C} 

V 

*i 

> + [K> 

V 

x i 

> = < 

V 

fi ► 


3b. 


3b. 


3b. 


f b. 


( 122 ) 


The RF transformation matrix [T] defined in equation (29) is altered slightly for this method of partition- 
ing: 
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(123) 


and is used to convert the above equation into new coordinates { q } p 
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, where 


[M],=[T] t [M][T], [C],=[T] t [C][T], and [K], =[T] T [K][T] . 


(124) 
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The transformed mass and stiffness matrices will be identical to those obtained for the free response 
shown in equations (31) and (32). Rather than carry out the multiplication of the righthand side of 
equation (124) to determine an expression, the transformation matrix and the original load vector are 
both carried throughout the analysis until the generalized force is needed, at which time the matrix 
multiplication is carried out. To actually perform the forced-response analysis, the coupled system above 
is transformed to an uncoupled system using a set of truncated mass normalized system eigenvectors 
[<E>]j resulting from a free-response solution of equation (124). 

{q}i = [ <J) ]i{q } 2 • 025 ) 


The uncoupled system is in the form 

[!]{q} 2 + [ C 1 2{q} 2 + [a] i {q} 2 - [^>] i T [x] T {f } , (126) 

where [A] 1 is the diagonal matrix of the eigenvalues of the transformed system. The standard (viscous) 
modal damping assumption is used for [C 2 ]: 


[C] 2 =2.0(0[A]i- (127) 

In this analysis, £is assumed to be constant for all the modes, and the input force is phased with a struc- 
tural mode to simulate bladed-disk excitation. There are standard techniques available for obtaining the 
frequency-response solution for this set of uncoupled single dof systems. Because of the series of linear 
transformations performed on the solution, a complex quantity for the solution is desired rather than one 
in terms of magnitude and phase, as is commonly derived in texts on the subject and which cannot be 
easily transformed. 

To this end, first express the element of the original excitation vector {ff acting on the n’th dof as 
follows: 


fn = Pn e '^ Qt 3n ^ = Pn elQte 13,1 = p n e lQt (cosa n - i sma n ) = F n e l£2t > (128) 

where p n is the amplitude, Q is the frequency of excitation, a n is the relative phase of the excitation 
force, and F n is a complex number defined as shown. The righthand side of equation (126) is defined to 
be the generalized force { f gen } , which can be calculated directly as 

{fgen} = {Fgen>e iQt = WW {F}e iQt . (129) 

It should be noted that this step now requires more information for the RF transformation matrix [T] 
than the free response does. As seen in equation (123), the columns of the transposed matrix [T] t corre- 


70 



sponding to the new forces applied at the internal dof’s f| are needed for the matrix multiplication to 
calculate the generalized force. In addition, for final recovery of displacement locations not at the 
boundaries, the same information will be needed. Therefore, additional elements of the modal matrix 
[O] 1 and of the residual flexibility matrix [G ib ] must be saved from the eigensolution and residual 
flexibility calculation and used to create a partial [T] matrix using equation (123). Only the partition of 
[T] corresponding to the internal dof’s where load is applied, internal dof’s where final displacements 
are recovered (if different), and the boundary partitions are actually needed for the calculation of {f gen } 
using equation (129). 

Now, it is assumed that the solution q 2 has a harmonic solution of the form 


q 2 j=Vje 


i(Qt-'Pj) _ Vj e i ^ t p - i ' I 'j 


(130) 


where j=l,...,p is the number of modes retained in [<f>] j and 4j is the output phase of the j’th generalized 
dof. Substitution of equation (130) into (125) yields the following expression for the j’th generalized 
coordinate: 


-VjQ 2 e 
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e 1 ^ 1 + /LjUje i ' P ’ 
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= F e 
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(131) 


or 


(-12 2 + 2£©jQ i + Aj)i>e = F gen je iQt , 

where tOj is the j’th natural frequency (square root of the eigenvalue Aj). Solving for o yields 


t>j=- 




F 

r gen j 


r _ \ 2 

^ x 1 

12 \ 

+ 

2C — 

i 

v®jj 


l 




(132) 


(133) 


Now, substituting (133) into (130), applying the definition of Fj as described in (128) and f gen as de- 
scribed in (129), and canceling out terms in the exponent yields the solution for the generalized coordi- 
nate q 2 : 


q 2 = 
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(134) 
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The complex vector {q} 2 can now be transformed back to the original coordinates {x} using equations 
(29) and (124): 


{ x } = [T][<J)] 1 {q} 2 . (135) 

The absolute values of this complex vector are then calculated to obtain the physical displacement 
vector. 


This solution procedure was applied to the deterministic two-bladed disk as a test case. A system 
of four harmonic loads, each of amplitude 1 lb and applied to be in-phase with a peak displacement 
location of mode 3, shown in figure 26, was considered. A range of excitation frequencies ± 5 Hz from 
the mode 3 deterministic natural frequency of 455.35 Hz was applied, and an equivalent viscous damp- 
ing value £ of 2 percent was used in all modes. The maximum resulting response was 0.044558 in. for 
node 640 (at the tip of blade A), occurring at the natural frequency, as expected. The deterministic 
solution procedure was verified by performing the same analysis in NASTRAN for the full-up model, 
using a frequency range of ± 5 Hz from the NASTRAN mode 3 frequency of 456.035 Hz; this run 
resulted in a value of 0.044832 in. for the same node, an error of only 0.6 percent. 



Figure 26. Two-bladed disk, mode 3 at 456 Hz. 
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It should be mentioned that an alternative, more direct method to solve for the vector { q } j could 
also have been applied. In this method, the matrices [M] t , [C]j, and [K]j, are not transformed into 
modal coordinates but are rather combined in a system matrix [Z] 1 as shown below: 


where 




(136) 


[Z], = -Q 2 [M], 


+ iQ[C] 1 +[K] 1 


(137) 


and f(to) is a complex quantity representing the magnitude and phase of excitation as shown in the 
derivation of the modal superposition method. The response vector { q(£2) } l is obtained by 
premultiplying f(£2) by the inverse of [Z] i at each excitation frequency of interest. This approach is 
preferred if there is no modal truncation performed, which is the case for the solution of {q} 2 after 
equation ( 1 26) is formulated. Since it is more direct, this method is recommended for future research. 

5.7 Forced Response — Probabilistic Case Using PDS 

As discussed above, applying the forced-response solution requires more information from the 
modal data and the residual flexibility matrices. This increases the number of dynamic rv’s necessary to 
solve the problem, which increases the complexity of the problem in several aspects. More statistical 
calculations, a larger correlation matrix, and more dimensions for the solution of the AFORM algorithm 
are all required. To minimize this number, it is important to decide a-priori which internal dof’s will 
either have external load applied or require a displacement recovery. Although the entire modal matrix 
and residual flexibility matrix are calculated, only those partitions of the modal matrix and the residual 
flexibility matrix, along with the boundary dof’s, are stored. In addition, only the correlations with these 
additional dof’s are generated. Note that one might also include the measured modal damping ratios of 
each substructure in the set of dynamic rv’s for the forced-response problem. Given the variability in 
measured damping properties, the incorporation of statistics on the modal damping value could be an 
important avenue for future research. 

To address the design problem defined earlier, the response variable for the forced-response 
analysis was chosen to be the maximum response for all dof’s in the structure. Initially, the excitation 
was applied at only a single frequency, but this was changed to a wide excitation bandwidth since the 
maximum response occurs at the damped natural frequency of the mode shape being excited, and this 
mode shape can occur over a range of frequencies for a probabilistic structure. It is assumed that the 
actual excitation mechanism could also vary in frequency by this amount, which is generally the case in 
engine turbomachinery, for instance. Using the maximum over some frequency range also helps to 
reduce the extreme variability in response which can be observed at any particular excitation fre- 
quency. 54 This variability, which is largest near peaks and zeros in frequency response functions, would 
introduce substantial nonlinearities in the response limit state surface. A range of ± 10 Hz about the 
deterministic natural frequency of 458 Hz was used in this case. This value is probably inadequate, since 
the free-response case shows that the actual range is between 438 and 472 Hz, so this may be a source of 
error in the final results. 
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The MC baseline frequency-response solution required extensive alteration of the MSCPOST 
subroutine within the NESSUS code. This subroutine reads the NASTRAN dataset OUGV1, which 
contains any displacement data obtained during the NASTRAN run. The frequency response solution for 
the baseline full-up model was run using the NESSUS/NASTRAN MC solution sequence, as with 
previous cases, but the entire displacement set for each sample was then read and a sorting routine 
performed to find the maximum value. The value was then fed back into the NESSUS main program for 
statistical calculations. This MC baseline solution sequence was extremely computer intensive, so only a 
limited sample size was run. 

The PDS procedure using both the MC and reliability approaches for the frequency-response 
problem is similar to the free-response problem. The main additions are the subroutines necessary for 
calculation of the different [T] matrices for each substructure and the complex solution algorithm de- 
tailed in the deterministic case. Because the value of the largest responding dof is no longer at a fre- 
quency known a-priori (since the third natural frequency and third natural mode will vary for each 
statistical sample), a sorting algorithm to find the maximum value was used that scanned the response 
values for all the selected dof ’s for all the frequencies in the bandwidth chosen. The dof ’s selected were 
those at the blade tips, which were assumed to contain the maximum responding dof for the entire 
structure with the chosen excitation. This value was then used as the NESSUS response variable, and 
input along with the corresponding load set. Also, to magnify the response a bit from the deterministic 
case, a modal damping ratio value of 0.5 percent was used. PDS MC and FORM cases were run as well 
as a normal distribution using the linear “approximate statistics.” The PDS quadratic AFORM method 
using the medians was also attempted, but the AFORM algorithm in NESSUS/FPI never converged to a 
solution. Several different attempts were made using different sizes of the variation of the input rv, but 
answers were not generated. The reason for this problem is unknown at this point, and is an issue for 
future research. 

The CDF results for the PDS methods are plotted along with a 1 ,000-sample MC baseline simu- 
lation in figure 27. The MC baseline CDF is extremely unsymmetric, so this clearly causes a problem for 
the reliability methods. It is also noted that there is fair agreement between the MC baseline result and 
that obtained from MC using dynamic rv’s. The errors for the design point at 99 percent are still around 
10 percent for both PDS methods, which is within reason especially when one considers that at present, 
there are no alternative methods for obtaining this maximum value. As with the free-response analysis, 
the significant error in the standard deviation and the lack of coincidence of the curves indicated visually 
signifies that there is still error present in the methods. An additional source of error in the reliability 
methods for this analysis is the extent to which the maximum responding dof does not have a smooth 
dependence on the rv’s. Possible discontinuities in this dependence would cause errors in the response 
surface generation created using numerical differentiation, as discussed in chapter 3. 

An examination of the computer run times shows that the MC baseline took 728 sec of CPU 
time and 27.3 hr of wallclock (which was much faster than normal due to low computer usage during the 
entire run), the PDS MC took 389 CPU-sec and 1,112 wallclock sec (the free- and forced-response were 
actually performed simultaneously), and the PDS FORM took 20.35 CPU sec, 60.6 wallclock sec, and 
less than 30 sec for the NESSUS/FPI run. Because these results showed an acceptable level of error 
from the design perspective and because the main goals of the method were achieved, an analysis of a 
realistic bladed-disk system was pursued and is discussed in the next chapter. 
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Figure 27. Two-bladed disk frequency response; MC versus PDS. 


6. FREE- AND FORCED-RESPONSE ANALYSIS OF REALISTIC BLADED-DISK SYSTEM 


6.1 Introduction and Literature Survey 

The development of the PDS method discussed in the previous chapters is now complete enough 
for it to be applied to the dynamic analysis of rocket engine or gas turbine disks with turbine blades 
attached. Because of finite manufacturing capabilities, which cause the geometry of the blades to be 
uncertain, and because of variability in the material properties, each of the blades typically have slightly 
different modal characteristics when considered independently; each of the dynamic characteristics can 
therefore be characterized by a distribution type, mean, standard deviation, and other statistical param- 
eters. An example of this variation is the first three measured cantilever modes of the turbine blades used 
on a first stage disk in the new Pratt & Whitney Space Shuttle Main Engine (SSME) high pressure fuel 
turbopump, shown in figure 28. The cov’s of these distributions are around 1 percent, and the same 
distribution fitting program as used in the benchmark problem was run on the first mode, yielding a close 
fit with either the lognormal or normal distributions. The blade-to-blade variations in dynamic character- 
istics result in the combined bladed-disk system being mistuned with respect to an ideal system having 
identical, evenly-spaced blades. It is known that the dynamic response of the structure becomes irregular 
and difficult to predict. Two of the main results of this irregularity are localization, which is the concen- 
tration of the response at a certain blade within the set; and amplification, which is the magnification of 
the response of the maximum-responding blade in a mistuned bladed disk over the maximum-responding 
blade (or blades) in a tuned disk. 

The mistuning phenomena has been examined in detail by researchers since D.J. Ewins described 
the problem in 1968. 55 Generally, two approaches have been taken for examining the problem. The first 
one is the deterministic approach, which uses both analytical and numerical techniques for deterministi- 
cally defined mistuned disks. A general review of the research by Leissa in 1981 summarizes some of the 
analytical methods used, such as the receptance method, in which elements of an impedance matrix 
relating force applied to the disk to response at the blade and vice versa are determined by using closed- 
form partial differential equations. 56 Ewins applied this method to the study of packeting of blades on a 
disc, where blades are split into groups with identical blade frequencies within the group. 57 MacBain used 
impedance to determine the approximate level of response for any particular level of mistuning, or level 
of spread of blade natural frequencies about the mean value, for a set of blades 58 Afolabi used these 
methods to define the localization and amplification phenomena seen in mistuned bladed disks. 59 

One of the major deterministic methods of examining these systems is to use a lumped mass 
model composed of N single-degree-of-freedom (sdof) spring-mass systems each grounded but connected 
to each other with a coupling spring to represent the bladed-disk system. These lumped masses can be 
used to generate sets of simultaneous equations for matrix solution. One of the first papers to identify 
localization in structures was published by Pierre and Dowell in 1987, who generalized a geometrically 
periodic structure with lumped masses; they also incorporated perturbation techniques to analytically 
model the mistuning level rather than use a given set of mistuned frequencies. 60 Pierre and Wei continued 
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Figure 28. Variation of first three first-stage turbine blade frequencies for P&W SSME fuel pump 






this localization work and examined the relationship with structural coupling. 61 They used a determinis- 
tic, closed-form solution for the eigenvalues and eigenvectors as a function of a single coupling 
parameter R and the “mistuning strength” of each sdof “blade,” which was defined simply as the varia- 
tion of that blade’s natural frequency from the mean, Af ; . The deterministic cases examined were 
generated by using a random number generator to determine the N Afj values. Only two complete 
mistuned systems were able to be examined for amplification and localization using this method. In 
addition to the lumped mass method, some finite element studies have been performed to examine 
deterministic mistuned bladed disks. Irretier created a simple deterministic finite element model of a 
mistuned bladed disk, similar to that described in the last chapter, and proved that CMS methods could 
be used to accurately reduce the size of the model and maintain the localization behavior of the sys- 
tem. 62 


The second major approach to examining mistuning phenomena is to apply statistical techniques. 
As with the deterministic approach, both analytical and numerical methods have been used. Generally, 
the perturbation method, as described in chapter 3, has been used to obtain the statistical properties of 
the response variable in the analytical approach. A closed-form expression for the variation of the ele- 
ments of the system matrix A, as expressed in equation (88), in terms of the input rv’s is obtained di- 
rectly. This is possible if the variation in the system matrix is kept extremely simple, as is the case when 
the system is composed of sdof “blades” and the only rv is the stiffness of each “blade.” 

One of the first statistical examinations was performed by Huang in 1982, who approximated 
the system as a continuous disk with statistically characterized randomly varying circumferential param- 
eters 63 A partial differential equation was developed for this continuous system, which was then solved 
by applying the perturbation method. The result was analytical expressions for the stochastic dynamic 
characteristics of the structure. Sinha and Chen built upon this work in 1989 by applying it to discretely 
modeled bladed disks using the sdof per blade approximation. 64 The perturbation method was used to 
create an analytical expression for the maximum responding blade of the structure, which could then be 
represented by a power series. The statistics of this series were then used to obtain probability density 
functions of the response. The expressions obtained in this paper are extremely confusing, and direct 
integration of the joint pdf is required to obtain the CDF, which limits the applicability to simple func- 
tions with a small number of rv’s, as discussed in chapter 3. Purely numerical statistical work has also 
been performed on the problem, notably by Griffin and Hoosac in 1983, who characterized the response 
of large lumped-mass systems by using an MC simulation to run multiple eigensolutions rather than 
determining an explicit expression for the response variable. 65 

Towards the beginning of the 1990’s, it became better recognized that a statistical approach was 
required for mistuning analysis because the blade population is large enough to require its characteriza- 
tion as an rv. In 1989, Pierre and Wei pointed out that the mean and standard deviation of the largest 
amplitude response experienced by a bladed disk cannot be generated by analytical techniques. 66 In this 
paper, they combined analytical techniques with MC simulations to understand the effects of mistuning. 
In a paper the next year, Pierre notes that only two studies of localization using statistical methods had 
been performed to that date. 67 In the same year, Hamade and Nickolaidis presented the first application 
of the reliability method of probabilistic analysis to the localization problem. They applied the FORM 
and the Rackwitz-Fiessler most probable point search algorithm to obtain the statistics on a dynamics of 
a simple two-span beam problem. Chamis and Shah in 1990 applied NESSUS to the bladed-disk prob- 
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lem, modeling the system with simple lumped masses with Young’s Modulus as the single rv. This study 
was just a cursory look at the problem as an application of NESSUS. In 1993 Mignolet and Lin ex- 
panded upon the closed-form expressions developed by Sinha, and made comparisons between mistuned 
analysis methods. 68 The comparison of the various curves were performed on a qualitative, visual basis. 

In 1995, Pierre and Kruse starting working on methods that would examine the statistics of more 
realistic mistuned bladed disks, realizing that the lumped-mass parametric studies previously performed 
are “difficult to relate to more descriptive finite element models.” 69 They emphasized that the use of 
MC simulations is “critical” in calculating the maximum response of bladed disks although it is “pro- 
hibitively” expensive. Therefore, along with Ottarson in 1995, 70 they developed a “reduced-order 
method” (ROM), a finite element reduction method similar to CMS that is aimed directly at bladed disks 
by modeling only one sector of the disk. Only the Young’s Modulus of each blade is set as an rv in this 
method, thereby allowing the natural frequencies of the blades to be directly mistuned. A deterministic 
tuned bladed-disk test case resulted in errors of 7 percent and 13 percent for the maximum blade re- 
sponse and system natural frequencies, respectively, between the “reduced-order” model and the full-up 
model where the size of the reduced-order model was 1/40 of that for the full-up model. A 1,000-sample 
MC simulation using a mistuned reduced-order model was performed to generate plots of the mean 
mistuned response normalized to the tuned response and a mean +3-0 mistuned response normalized in 
the same way. Using this information, the researchers were able to identify important characteristics for 
mistuned bladed disks, in particular that “moderately weak” blade coupling is required for significant 
amplification. Yang, Griffin, and Kiefling simultaneously developed a very similar reduced-order ap- 
proach, but did not perform any statistical analyses. 71 

6.2 Comparison of Previous Techniques for Mistuning Analysis With PDS 

The extensive literature survey summarized above clearly identifies the importance and unique- 
ness of the PDS technique, not only for probabilistic substructures in general, but in particular to the 
analysis of mistuned bladed disks. One of the biggest drawbacks to most of the studies is that they 
concentrated on very simple systems. The SSME turbine blade, for instance, is a doubly curved, hollow 
structure with 25 modes in the excitation range; a typical finite element model, shown in figure 29, 
contains over 30,000 dof ’s. Representing this structure with an sdof spring-mass system clearly is 
inadequate for most analyses. For this reason, while the analytical, closed-form solutions as developed 
by Sinha and Mignolet are important for understanding some of the physical phenomena of the problem, 
they are impractical for realistic analysis. At this point in time, practical analysis is vitally needed by 
industry. Rocketdyne, the manufacturer of the standard set of turbopumps for the SSME, performed a 
quick analysis of their second stage liquid oxygen pump and showed that considerable amplification did 
exist for a system model as compared with a single blade approach. 72 They reported in 1989 that they 
and other bladed-disk manufacturers are not using any kind of mistuning in their production analysis 
because the research in the field has been too unrealistic. In addition to this practical need, the analysis 
of realistic structures is important to capture the true characteristics of the mistuning phenomena, which 
can be lost using simplified models. 

The work performed recently by Pierre, Griffin, and their co-workers does attempt to move 
towards realistic analysis methods. However, because of their reliance on the perturbation method, 
which is only practical for blades represented by sdof structures where the partial differential of the 
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elements of the system matrix [A] with respect to each rv can be represented by closed-form equations 
(as discussed in section 3.5), these investigators have not taken advantage of the tremendous advances 
made in probabilistic analysis using reliability methods. They therefore still model the variations in the 
blades by a single rv, which is clearly overly restrictive. Sinha briefly referenced the FORM in the mid- 
1980’s, but quickly dismissed it and has not examined it since. Hamade and Nickolaidis used the method 
for localization, but have not done any work since then, extending it to bladed disks; and Chamis 
touched on the subject, but has not performed any additional work. 



Figure 29. Finite element model of SSME turbine blade. 

The PDS method has a number of advantages over the methods discussed in the preceding 
survey. First, realistic finite element models can be used for two reasons: the size of the model can be 
reduced tremendously using component mode synthesis, and each blade on the disk, which can number 
up to a 100, does not have to be modeled separately since the statistical characteristics of each blade are 
identical. A second advantage is that assigning the dynamic characteristics of the structures to be the rv’s 
makes no assumptions regarding the origin of the random characteristics of the blades; all of the preced- 
ing research assumes all of the randomness is in the natural frequency or elastic modulus of the blade, 
and completely ignores variability in the mass, geometry, etc. which may only show up in the mode 
shapes or residual flexibility. Using a realistic model also allows actual values of mistuning and coupling 
to be inherent in the system model rather than having to specify a “mistuning” or “coupling” parameter. 
Finally, applying the approximate reliability techniques to obtain the response statistical density func- 
tions is an order of magnitude faster than MC, which is imperative for the completion of parametric 
studies necessary for understanding these types of complex systems. 
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6.3 Application of PDS for Nine-Bladed Disk 


6.3.1 Model Definition, Baseline Analysis 


The application of PDS on a realistic nine-bladed-disk system, shown in figure 30, can now be 
attempted. A great deal of parametric investigation is necessary for a complete mistuning analysis; the 
goal here is to obtain some basic results that illustrate the applicability of PDS to the problem, thereby 
laying the groundwork and methodology for future studies of realistic systems. The disk and blade “A” 
of the model are identical to the system described in chapter 5, and the other eight blades possess identi- 
cal yet independent statistical parameters to blade “A.” Therefore, the construction of only one blade 
model is required and rotation transformations can be applied to this model to create the bladed-disk 
system model. Compared to some other methods of modeling mistuned bladed disks in which each 
slightly different blade is modeled separately, these rotations are a tremendous savings in effort and 
computer time since the substructure has already been reduced using CMS. 
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Figure 30. Nine-bladed disk, mode 6 at 444.6 Hz (using CMS). 
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The first step in the analysis is to produce baseline results for comparison. This was accom- 
plished using the same free- and forced-response procedures as described in chapter 5. Based upon the 
close agreement in CDF’s between the 400- and 1,000-sample cases for the two-bladed-disk system, 400 
samples were used for this larger system to save computer time. The assumption that this number of 
samples is sufficient may not be valid and could be a source of error. As discussed in the previous 
chapter, this will be the first publication of an MC analysis of a realistic bladed-disk system modeled to 
this degree of detail in both modeling representation and level of stochasticity. 

Following the PDS procedure, an MC simulation using NESSUS/ NASTRAN to obtain the 
statistics of the first blade was then performed. The mean and standard deviation of the resulting distri- 
butions for the dynamic rv’s and the correlation matrix were then used to create independent “load 
cases” used for creating the substructure stiffness matrices for coupling with the system. Only those 
uncorrelated rv’s whose correlation matrix eigenvalues were greater than 3 percent of the sum of the all 
values were actually used as dynamic rv’s in NESSUS, as previously discussed. This approach was 
taken to eliminate eigenvalues which did not contribute to the total randomness in the system, as repre- 
sented by the sum. Rather than creating an enormous load case matrix [u] containing independent 
variations of each rv for each blade, a matrix [u] for the first blade only was created and reused for all 
the blades. This was accomplished by coupling each variation case for blade i=l,...,9 with the mean 
value case (a column of all zeroes) for all the other blades. A rigid body transformation matrix 


cos0, sintL 


R = 


'z 

-sin# 

0 


COS# 


0 

0 


z 

0 1 


(138) 


was created, where 0 Z is the angle about the z axis of the original blade from the transformed blade. In 
degrees, the blade angles are: 

# z =0, -40, -80, -120, -160, -200, -240, -280, -320 . (139) 

These rotation angles were then used on each node in the eigenvectors and the residual flexibility matrix 
by placing [R] into partitions along a diagonal band in a global transformation matrix. The eigenvalues 
are the same for each blade, and the correlation matrix is the same also because only the z translational 
displacement elements of the above matrices are used as dynamic rv’s, and the rotation about the z axis 
will not affect these rv’s. The substructure matrices resulting from applying the RF CMS method using 
this information was then coupled into the system using a map identifying the location in the system 
matrices of the boundary dof’s, the dof’s whose displacements were chosen to be recovered, and the 
generalized dof’s. 

6.3.2 Free-Response Analysis 

An analysis to obtain the CDF of mode 6 (see fig. 30) was performed on this nine-bladed disk. 
This mode was chosen because it was used as the excitation mode for the frequency response analysis 
performed subsequently. Both PDS MC and PDS reliability methods were applied. As described in the 
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last chapter, the medians of the dynamic rv’s were substituted for the means in the variable transforma- 
tions. The results are shown in figure 31, with the figure indicating the degree to which the CDF curves 
match the baseline and the table quantitatively expressing the error using the various measures described 
in chapter 5. In general, the AMV update of the linear - FORM method yielded the best results, with a 
skew at the mean of only 1.2 percent, and errors at the CDF values of interest for design (0.01 and 0.99) 
of 2.5 and 0.6 percent, respectively. The shape of the CDF does show some error as indicated visually 
and by the standard deviation error of 48 percent; this average standard deviation value for comparative 
purposes was calculated by taking the difference in the response values at ±lc?(0.8414 and 0.1586) and 
dividing by 2. The PDS MC results actually show worse results in both the shape of the CDF and in the 
error measurements. This result was unexpected since, as described in the previous chapter, there are 
fewer sources of error in the MC methodology than the PDS reliability methodology. For both the free- 
and forced-response analyses, it should be noted that statistical outliers (values beyond ±4a) generated 
in the PDS MC analyses were deleted from the sample to not falsely exaggerate the standard deviation. 
The APORM method yields the worst results. A possible reason for these results is discussed in the next 
section. However, the approximate statistics generated by the partial quadratic formulation of the limit 
state in the AFORM solution was used to create a CDF by assuming a lognormal distribution. This CDF 
actually matches the AMV results almost exactly. These results are consistent with a conclusion that the 
function is nonlinear without having a large coefficient of variation. 
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Figure 31. Nine-bladed disk, mode 6 PDS evaluation. 
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Examining the run times shows the MC baseline using 400 samples took 1,527 sec of CPU time 
and 23.2 hr wallclock, the PDS MC using 500 samples took 3,867 CPU-sec and 4.8 hr wallclock, and the 
PDS FORM took 778 CPU sec, 2,849 sec (0.79 hr) wallclock, and < 30 sec for the NESSUS/FPI run. As 
expected, the PDS methods result in a tremendous savings in overall computer time. 

There are several interesting points to note in the results of the free-response analysis that could 
explain some of the errors. As with the two-bladed-disk system, a deterministic NASTRAN analysis of the 
system was run using the means of the primitive rv’s as well as a deterministic RF CMS analysis using the 
medians of the dynamic rv’s. The NASTRAN deterministic analysis yielded values of 445.6 Hz for both 
modes 6 and 7, and examination of the mode shapes verified that their only difference is a rotation about 
the center of the disk. However, the median of the MC baseline run for mode 6 was 439.6 Hz, a decrease 
of 6 Hz. For a tuned system, the NASTRAN deterministic result would be almost identical to the MC 
median value (for enough samples), but because of the mistuning, the double mode at 445.6 hz bifurcates 
into a low mode (number 6) and a high mode (number 7), and since the lower mode is being sampled, the 
CDF will be skewed lower. The RF CMS deterministic analysis using the dynamic rv’s, on the other hand, 
results in mode 6 at 444.87 Hz and mode 7 at 451.2 Hz. Apparently, some aspect of the reductions in the 
RF CMS methodology, which would be expected to stiffen the system, actually slightly mistuned the 
system, causing the double mode to split for this median case. The median (CDF = 0.5) value for the 
FORM analysis is also defined to be this value, so the skew between the FORM CDF and the NASTRAN 
MC baseline CDF at the median can be explained by the reduction of the baseline due to the bifurcation. 
The increase in the median value for the PDS MC results compared to the RF CMS deterministic value, 
on the other hand, cannot be explained using this reasoning and is not understood at this time. 

6.3.3 Forced-Response Analysis 

A frequency-response solution for the maximum responding blade tip dof in response to an excita- 
tion wave in the three nodal diameter shape of mode 6 and 7 at a deterministic mode 6 frequency of 444.6 
Hz was then run using the synthesized systems. This solution yields a more general evaluation of the PDS 
technique than the analysis of mode 6 only because, as described in chapter 5, a modal superposition 
method is used so the effect of the variation in all the modes is included, although the most important ones 
are modes 6 and 7. A 2-percent modal damping value was used in the analysis. After several iterations, a 
frequency band of 50 Hz was determined to be required for the accurate determination of the peak re- 
sponse since the mistuned system can resonate at both modes 6 and 7, which are split by the mistuning. 
The maximum value of the response was obtained by sorting through the responses of all the blade tips at 
all analyzed frequencies. For every “load case,” the dof and frequency of the maximum response was 
output. As expected, this spatial and spectral location varies significantly, both due to the variation in the 
mode 6 and 7 frequencies and the mistuning localization effects. 

The reliability methods applied for the free-response analysis were applied in this case as well. In 
addition, the AMV+ method, which is a further update method in which the limit state is expanded about 
each most probable point, was used. However, none of the higher order methods achieved reasonable 
results, and so are not included in the summary chart, which is shown in figure 32. Note that there is good 
agreement between the shape of the MC baseline and MC PDS results, although the median values differ 
by about 5 percent. The fact that the curve for the PDS MC is shifted to lower levels of response can 
partially be explained by the stiffening that accompanied the RF CMS model reduction, which was also 
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seen by the higher frequencies of modes 6 and 7 in the previous section. The relatively good agreement 
in median values predicted by the PDS MC and PDS FORM methods is also reassuring since both 
utilize the same reduced models. The response levels at a CDF value of 0.99, which is of significant 
importance in the design process, show that the FORM has an error of 6.25 percent and the PDS MC an 
error of 14.13 percent. Examining the run times shows the MC baseline using 400 samples took 780 sec 
of CPU time and 23.0 hr wallclock, the PDS MC using 500 samples took 14,627 CPU sec and 4.16 hr 
wallclock (the wallclock time for this run was unusually close to the CPU time because it took place 
over a period of very low computer usage), and the PDS FORM took 3,088 CPU sec, 2. 19 wallclock hr, 
and < 30 sec for the NESSUS/FPI run. These results indicate that the PDS methodologies can provide 
the designer with a legitimate tool for generating previously unobtainable maximum response values for 
a mistuned bladed disk in a reasonable amount of time. 

To determine whether the mistuned system analyzed actually does display amplification effects 
as described in the bladed-disk literature, an analysis of a tuned system was also performed. To obtain a 
true “apples-to-apples” comparison, this analysis must also be probabilistic, but the values of the rv’s 
within each blade must be identical to those values in the other blades for the disk to be tuned. An MC 
baseline run was therefore performed on a full-up model with only three independent rv’s for the whole 
system, which are the same rv’s as for the mistuned system but set to be the same for all the blades. As 
alluded to previously, the deterministic response of the system using the median values of all the input 
rv’s would be for a tuned, (not mistuned) system; a median mistuned system is difficult to quantify other 
than with a probabilistic analysis. Figure 33 compares the tuned and mistuned MC CDF plots. Also 
shown is the response amplitude for the deterministic case (small circle). This plot clearly shows an 
amplification of the response, indicating that the model successfully represents a mistuned system. 

6.3.4 Discussion of Error and Conclusions 

A possible reason for the problems with the approximate high order methods for both the free- 
and forced-response solutions is the extreme nonlinearity of the response away from the mean value. 

This nonlinearity is determined by isolating the response as a function of the two independent rv’s that 
NESSUS indicated were the most important to the response under its “sensitivity” column, which is part 
of its standard output. The response to these two variables, which turned out to be numbers 23 and 89 
out of 99 total independent rv’s for the system (the first rv’s for the third and ninth blade), was calculated 
for a range of -2 times the standard deviation to +2. Plots of the value of the mode 6 natural frequency 
and of the maximum blade response versus these variables is shown in figure 34. The maximum blade 
response is jagged, which may be due to its possible discontinuities as a function of the input rv’s, as 
discussed in chapter 5. It is apparent that the function cannot be approximated by a quadratic; a linear 
approximation about the mean may actually be a better fit of the response surface. For the mode 6 
frequency response, the minimum inflection point appears to be somewhat less than the mean, so a 
quadratic approximation about the mean would also not be very accurate. It must be remembered that 
the independent rv’s are not tied to a specific primitive rv; therefore, it should not be surprising that the 
curve is nonmonotonic. It may be noted that expanding the limit state about a value equal to -0.3 times 
the standard deviation might give good results for a quadratic approximation technique. It is recom- 
mended that studies along this line be pursued in future research. 
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Figure 32. Maximum frequency response using PDS and MC. 



Figure 33. Frequency response of maximum blade to 3ND input for 50 Hz frequency range at mode 
6 of median system; comparison of tuned system variation versus mistuned system. 
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Because of the unusual phenomena of closely spaced natural frequencies, localization, and 
amplification, it was not known if a reasonable result using the PDS methodology would be possible for 
the analysis of a mistuned bladed disk. The results shown in this chapter do show significant discrepan- 
cies with MC baseline runs, but they are promising enough to be a valuable tool for analyzing realistic 
bladed disks in new engine designs. This is especially true if one considers the accuracy, availability, and 
computational intensity of other analysis techniques. Further refinement of the method may be necessary 
to use it for exploration of the whole range of mistuning issues. 




Number of Standard Deviations of Independent rv 


Figure 34. Nonlinear dependence of mode 6 and maximum response as a 
function of independent rv’s 23 and 89. 
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7. CONCLUSIONS AND FUTURE RESEARCH 


7.1 Conclusions 

A methodology has been developed in this thesis to perform probabilistic analysis of structural 
systems composed of substructures whose dynamic characteristics can be characterized statistically. This 
method, entitled probabilistic dynamic synthesis, has been created using new developments in probabi- 
listic analysis and the residual flexibility method of component mode synthesis. It enables the use of 
realistic, comprehensive statistics obtained from modal testing of substructures rather than much less 
complete statistics for primitive rv’s such as geometry variation. These statistics can be applied either 
using the well-tested but computationally intensive method of MC simulation, or by using new reliabil- 
ity methods, which can significantly reduce the computational effort. Additionally, the use of component 
mode synthesis to couple these substructures provides for a considerable savings of computer resources 
through a reduction in the size of the analytical model of the structural system. 

The development of this methodology initially required the examination of the many methods 
available for component mode synthesis. Investigations into the final makeup of the system matrices 
resulted in the conclusion that the residual flexibility method was the only one that solely required 
information that could be obtained from test. This allowed the complete statistical variability of the 
substructure population to be accounted for in the subsequent analysis. 

An evolution of the major developments in probabilistic methods was then presented to identify 
and understand how best to apply these methods. The approximate “reliability” methods were found to 
be able to incorporate realistic structural analysis techniques such as the finite element method, and can 
provide accurate CDF’s using an order of magnitude less calculations than MC simulations. Methods 
using both first order and partial second order approximations of the limit state function were found to 
be applicable to the PDS method. In addition, a detailed comparison was made between the perturbation 
method, which is an analytical method used frequently in the study of periodic disordered systems, and 
the reliability methods. Essentially, the perturbation method was found to be equivalent to the FORM 
but much less applicable to realistically sized systems. 

A simple benchmark test problem consisting of two probabilistic substructures was then ana- 
lyzed. Many of the details involved in applying the PDS method to a structure were developed using this 
test problem, and a variety of reliability methods were applied to obtain statistics of the system funda- 
mental frequency. A procedure was developed that used MC numerical analysis to simulate the modal 
testing that would actually be used in PDS to determine the statistics of the substructure eigenvalues, 
eigenvectors, and residual flexibility matrix elements. These statistics were then used for the generation 
of “load” cases for a derived set of independent, un correlated rv’s that were used to obtain the desired 
response variable in terms of variation of each of the independent rv’s. These “datasets” defining the 
response surface were then input into a newly available probabilistic code NESSUS which uses reliabil- 
ity methods to determine the CDF of the response variable. For this simple case, a linear approximation 
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of the limit state using the first order reliability method provided results that matched well with a com- 
parison to an MC simulation of the full, unsubstructured system. For larger coefficients of variation 
(10-15 percent), the partial higher order advanced first order reliability method resulted in better an- 
swers. The advanced mean value method was also applied, and it was found that this method is better 
suited to problems where the limit state itself is more nonlinear. 

The PDS method was then applied to a realistic finite element model of a three-substructure 
system. The baseline structure was modeled with a sufficient number of independent primitive rv’s to 
accurately represent the random variation in realistic structures, which had not been done by other 
researchers. A great many practical issues necessary for actual use of the method were discovered and 
resolved, including developing the ability to perform an accurate baseline MC analysis of a realistic 
structural system, reducing the number of independent rv’s by examining the relative size of the eigen- 
values of the correlation matrix, and using the medians, rather than the means, of the samples of the 
dynamic rv set for the transformation to the independent rv set. These dynamic rv’s were applied using 
both the MC method and the linear and higher order reliability methods. The accuracy of the results was 
sutficient for design purposes and showed a tremendous savings in computer time, but there were 
discrepancies in the CDF’s that should be investigated further. 

In addition to the solution for a selected eigenvalue, the statistics of the forced, frequency- 
response solution for the system were obtained. This required the derivation of a methodology for 
applying phased excitation in the frequency domain and implementation to the residual flexibility 
method of CMS. The errors for the PDS methods were larger than the free-response case, and there were 
substantial differences in the shapes of the CDF’s from the baseline case. From a design viewpoint, 
however, the errors were within the acceptable range. Considering the lack of alternatives for probabilis- 
tic forced-response analysis and the tremendous savings in computer time, the PDS method was con- 
cluded to be an important new tool for analyzing these types of systems. 

The PDS method is particularly suited for the mistuned blade-disk problem because of the 
random nature of the substructure blades, so it was then applied to this problem. An extensive literature 
survey was performed that clearly identified the advantages of the method over other techniques used, 
such as the perturbation method. The motivation for analyzing realistic disks, both because of the practi- 
cal nature of the problem and to gain a better insight into mistuning, was also presented. A disk with 
nine statistically identical but independent blades was modeled, and the previously developed procedure 
applied to obtain the free- and forced-response solution. Although a high degree of nonlinearity was 
associated with the problem because of the mistuning phenomena, answers similar to those obtained for 
the simpler three-substructure system were obtained for this larger bladed disk . 

In conclusion, the PDS method was shown to be an important new method for the determination 
of the statistical response of nondeterministic structures. Acceptable results for a variety of structures 
have been shown, and the details of the procedure have been examined. In particular, application of PDS 
to the mistuning problem should prove extremely beneficial for both an increased understanding of the 
phenomena and for actual design of bladed disks by industry. 
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7.2 Directions for Future Research and Application 


A great many avenues of research lie ahead for the application of this research and for its contin- 
ued development. To help serve as a roadmap for this research, table 7 is presented to clarify the meth- 
odologies examined and developed in this thesis, and to compare their benefits and limitations. This 
table also can be used to identify the options available to an engineer when designing a bladed disk or 
other structure with a significant probabilistic component. 

Several of the techniques developed in this thesis are of immediate benefit for the analysis of 
general probabilistic structures and bladed disks in particular. An immediate application of these meth- 
ods will be for the bladed disk presently being developed by NASA/MSFC for the turbopump of the 
new X-33 reusable launch vehicle. This information will be extremely important for design of the 
blades, which have been assumed to be tuned until now, as is standard practice in industry, with various 
safety factors applied to account for large responses. The use of the techniques and computer codes 
developed in this thesis can quantify the statistical distribution of the actual response values, providing 
much more accurate information for engineering design and analysis. 

While the PDS method is being applied to the above engineering problem, several options should 
be explored to increase its accuracy. One major source of improvement to the PDS reliability methods 
may be obtained by more direct integration of the FORTRAN codes written for this development with 
the automated features of the NESSUS code. This may eliminate some programming error and would 
allow use of more advanced reliability techniques, such as the AMV+ method, convolution method, 

Latin Hypercube sampling, and the FFT method. 73 Other improvements may be possible by reducing the 
other sources of error discussed in this thesis, such as using a larger sample size for the MC analysis of 
the original primitive rv’s, expanding the limit state about more suitable points than the mean, and 
incorporation of more modes in the reduced model of the bladed disk to ensure the calculation of double 
modes. Finally, a significant source of error for the forced response problem can be addressed by exam- 
ining only an sdof, rather than the maximum responding dof. This should increase the smoothness of the 
limit state and the accuracy of the reliability methods. A greater understanding of these errors may be 
obtained by generating the PDF of the results in addition to the CDF. The PDF would highlight the 
shape of the distribution, which has been shown to be the main discrepancy between the MC baseline 
and PDS results. 

Performing the PDS analysis on an engineering system may require some variation of the tech- 
niques developed in this dissertation to account for the unavailability of the data. Since residual flexibil- 
ity measurements may be difficult to obtain in certain circumstances, a hybrid test/analytical PDS 
methodology could be developed which uses the analytical value of the residual flexibility and measured 
eigenvalues and eigenvectors. As mentioned in chapter 2, the measurements of these quantities them- 
selves are not without error. Present testing methods actually require a statistical sampling of a single 
structure to generate an output value. The fact that the PDS method requires many samples, however, 
may actually mitigate this error because it uses the statistics of the measurements directly in the problem 
formulation, unlike a deterministic problem. An evaluation of the robustness of the statistics taken for 
PDS is an important subject for future research. 


90 



Table 7. Comparison of methods for analysis of structures with probabilistic substructures. 


Method, Name, 
Users, Year 

Model Size 

Random Variable 
Description 

Modeling 

Methodology 

Probabilistic 

Methodology 

Benefits 

Limitations 

Questions 

“Lumped mass 
approach:” Sinha, 
Pierre, etc. 1980— 
1990 

One lumped mass 
per substructure 

One primitive 
rv/substructure 

Closed-form 

equations 

Perturbation — 
closed-form 
equations of [A] 
wrt. rv’s 

Quick 

Inaccurate for real 
structures since neither 
1 dot or 1 rv/ 
substructure are 
correct 


“Reduced order 
methods:” Pierre, 
Griffin, 1995 

Finite element 
model 

One primitive 
rv/substructure 

Reduced order 
methods— similar 
to CMS 

Monte Carlo 

Representative structure 

(1) Unrepresentative 
description of rv’s; 

(2) Slow, CPU-intensive 
probabilistic calculations 


“Full model 
Monte Carlo” (MC 
baseline): in this 
thesis; previously 
unpublished 

Finite element 
model 

Multiple primitive 
rv’s/substructure 

Full model 

NESSUS/ 
NASTRAN Monte 
Carlo 

(1) Representative structure, 
more representative rv’s 

(2) Accurate CDF for given 
rv’s 

(1) Rv’s difficult to 
characterize accurately 

(2) Slow, CPU 
intensive 


“Full model- 
reliability:” 
immediate future 
work 

Finite element 
model 

Multiple primitive 
rv’s/substructure 

Full model 

Reliability 

methods 

(1) Representative structure 

(2) More representative rv’s 
than above 

(3) Faster than “full-model MC” 

(1 ) Rv’s difficult to 
characterize accurately 

(2) Faster probabilistic 
method than “full- 
model MC” 

May be less 
accurate CDF for 
given rv’s 

“PDS — Monte 
Carlo:” in thesis 

Finite element 
model 

Complete set of 
dynamic rv’s from 
test/substructure 

Residual 
flexibility CMS 

Monte Carlo 
(FORTRAN) 

(1) Representative structure, 
most representative rv’s 

(2) Accurate CDF for given rv’s 

(3) Faster than “full-model MC” 
(more efficient modeling) 

Slow, CPU-intensive 

probabilistic 

calculations 

Results not as 
accurate as 
anticipated 

“PDS— Reliability:” 
in this thesis 

Finite element 
model 

Complete set of 
dynamic rv’s per 
substructure 

Residual 
flexibility CMS 

Reliability 

methods 

(1) Representative structure, 
most representative rv’s 

(2) Most efficient model 
representation and 
probabilistic methodology for 
speed 

Less accurate CDF 
than “full model— 
Monte Carlo” for 
given rv’s 

Further improvements 
in accuracy req’d 
before use in general 
bladed-disk studies 

“CMS only- 
reliability:” future 
work 

Finite element 
model 

Multiple primitive 
rv’s/substructure 

Craig-Bampton 
component mode 
synthesis 

NESSUS/ 

NASTRAN 

reliability 

methods 

(1) Representative structure, 
fairly representative rv’s 

(2) Most efficient model 
representation and probabilistic 
methodology for speed 

Rv’s difficult to 
characterize accurately 





Another assumption requiring future investigation is accurate representation of the boundary 
connections. In a traditional bladed disk, the blades are inserted into slots in the disk and then held in 
place by a combination of centrifugal force and axial stops. In this case, the interface of the blade with 
the disk is difficult to model accurately. A study using the boundary connection as an independent, 
analytical rv may prove important for representing these structures. This problem is not as important for 
blisks, in which the blades are integrally machined with the disk. In these cases, the free-free modes 
would have to be obtained by either cutting the blades from the disk or by clamping the disk tightly and 
using the techniques developed by Han, Chen, and others for converting cantilevered modal test results 
to free-free results. 

A similar extension of the PDS method might be in the determination of the statistics of the 
global damping value of a structural system as a function of the damping values of the substructures. 
Modal damping ratios are generally measured during modal testing, so statistics could be compiled and 
the ratio used as a dynamic rv. These values could be used directly in the substructure forced response 
equations. Estimates of the global damping value for a particular mode could be obtained from the 
measured frequency response functions using, for example, the half-power method. This global value 
could then be used as the desired output variable of interest and a CDF obtained. 

In addition to the actual application of the PDS method to an engineering problem, it would be 
illustrative to quantify the gain in accuracy in PDS enabled by the use of the dynamic rv’s. These rv’s 
allow a very general degree of randomness in the structure compared to other methods that can only use 
a single rv per substructure, like Young’s Modulus. This could be accomplished by comparing the results 
of the MC baseline analysis presented in this thesis with MC results obtained for an identical bladed 
disk, but with only a single rv represented. The inherent inaccuracy of primitive rv’s has not been recog- 
nized in the literature, and the use of dynamic rv’s to model more representative random structures is 
one of the major achievements of this research. 

Finally, the PDS method can be used to examine the mistuning phenomena in greater detail. To 
this end, the variety of new tools and techniques developed and tested during this research, as shown in 
table 7, can be applied. The first of these to be applied would be the “full-model MC,” “full-model 
reliability,” and “CMS only — reliability” methods, which use full-up or substructured analyses of purely 
analytical mistuned bladed disks modeled with primitive rv’s and which use either the MC or reliability 
methodologies to obtain the desired statistical response. These analytical results would be used for 
parametric studies and then compared with results from using the “PDS MC” and “PDS reliability” 
methods, which are more applicable for an actual bladed disk. The most recent work by Ottarson and 
Pierre 74 has identified several important relationships concerning the amount of amplification as a 
function of damping, coupling, and mistuning strength, and these results could be verified for the more 
complex structures now able to be analyzed with these new tools. These results could help conclusively 
define the effects of mistuning for realistic bladed disks, which will be of great importance to both the 
industrial and academic communities. 
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